Transcription

1 R EZR *1 1 R R [ ] R Rgui EZR (cross-over design) F Welch t Wilcoxon Wilcoxon (Kruskal-Wallis) Fligner-Killeen Dunnett *1 1

2 [email protected] 1 R R MS Windows Mac OS Linux OS *2 Windows Mac OS Linux tar ubuntu R R CRAN (The Comprehensive R Archive Network) CRAN *3 *4 *5 1.1 R [ ] Windows CRAN R *6 *7 SDI (separate windows) Macintosh Mac OS X CRAN R pkg Tcl/Tk tools *8 *2 R *3 *4 *5 *6 R win.exe *7 R 1 MDI Rcmdr/EZR *8 2

3 Linux Debian RedHat ubuntu CRAN CRAN R tar.gz./configure make make install 1.2 R Windows Linux Mac OS X \ Windows R Rgui R *9.Rprofile.RData > > R + Windows ESC Source source(" ") Windows / \\ * 10 R bin Windows 7/8 R R R 1.3 Rgui q() < X X <- c(1, 4, 6) function() meansd() *9 (S) R_USER.Renviron R_USER="c:/work" proxy Windows proxy --internet2 R LANGUAGE="en" R MDI SDI --sdi *10 \ 3

4 meansd <- function(x) { list(mean(x), sd(x)) } install.packages() CRAN Rcmdr install.packages("rcmdr", dep=true) dep=true dependency Rcmdr Rcmdr TRUE T T TRUE? t t.test?t.test <<- 1.4 EZR R John Fox Rcmdr R Commander R Rcmdr EZR * 11 EZR Rcmdr install.packages("rcmdrplugin.ezr", dep=true) Rcmdr library(rcmdr) R Commander GUI R Commander library(rcmdr) Rcmdr Commander() detach(package:rcmdr) Rcmdr library(rcmdr) R Commander GUI EZR Rcmdr RcmdrPlugin.EZR OK R (Y) EZR Rcmdr * 12 *11 * Rcmdr RcmdrPlugin.EZR EZR R Windows R EZR 4

5 : 3 3 triplicate duplicate triplicate triplicate % 5% R.A. Fisher t 2 t t (RA) 5

6 2.2 (randomization) Fleiss JL (1986) The design and analysis of clinical experiments (random number table) (random permutation table) 2 Welch t 3 (One-way ANOVA) 2.3 (cross-over design) 2 Hilman BC et al. Intracutaneous immune serum globulin therapy in allergic children., JAMA. 1969; 207(5): A 34 2 B A B * 13 PID Ord.LD LD DL LD L D T0.Alb T1.Alb *13 web 6

7 PID Ord.LD LD T0.Alb T1.Alb 1 LD L LD L LD L LD L LD L DL D DL D DL D DL D DL D LD D LD D LD D LD D LD D DL L DL L DL L DL L DL L ) T0.Alb T1.Alb (0 ) LD Ord.LD OK Univariate Type III Repeated-Measures ANOVA Assuming Sphericity SS num Df Error SS den Df F Pr(>F) (Intercept) e-15 *** Factor1.LD ** Factor2.Ord.LD Factor1.LD:Factor2.Ord.LD Time e-08 *** Factor1.LD:Time e-05 *** Factor2.Ord.LD:Time Factor1.LD:Factor2.Ord.LD:Time Signif. codes: 0 *** ** 0.01 * Pr(>F) Factor1.LD 5% Time 5% Factor1.LD:Time 5% Factor2.Ord.LD Pr(>F) 2 7

8 3 60 kg, 66 kg, 75 kg R mean(c(60,66,75)) ( )/3 Microsoft Excel 10 ID (cm) (kg) Microsoft Excel R Excel *.xls OpenOffice.org calc *.ods Microsoft Excel (F) (T) (*.txt) xls txt (S) OK Excel desample.txt R Dataset 1 8

9 Dataset <- read.delim("desample.txt") Rcmdr a Dataset OK Excel OK Windows R Rcmdr1.4-9 RODBC Excel from Excel, Access, or dbase dataset b Excel a EZR URL b EZR Excel Access dbase 4 survey EZR MASS MASS survey OK survey 237 Sex Wr.Hnd cm NW.Hnd cm W.Hnd Fold 3 Pulse Clap 3 Exer 3 Smoke 4 Height cm M.I cm m Age 9

10 4.1 X R barplot(table(x)) EZR Smoke R fx <- table(x) barplot(matrix(fx,nrow(fx)),beside=f) Rcmdr EZR 100% R px <- table(x)/nrow(x) barplot(matrix(pc,nrow(pc)),horiz=t,beside=f) Rcmdr EZR 4.2 R hist() Sturges right=false R (Age) hist(survey$age) hist(survey$age, breaks=1:8*10, right=false) Rcmdr EZR survey Age Rcmdr EZR R qqnorm() survey (Pulse) qqnorm(survey$pulse) EZR QQ (Kolmogorov-Smirnov ) (stem and leaf plot) 5 10 R stem() stem(survey$pulse) 10

11 EZR Pulse OK Output (box and whisker plot) 1.5 R boxplot() survey (Smoke) (Pulse) boxplot(survey$pulse ~ survey$smoke) EZR survey 0 1 Smoke 1-1.5x x OK survey Smoke Pulse (scatter plot) R plot() plot() pch points() symbols() matplot() matpoints() pairs() text() identify() survey R (Age) (Height) plot(height ~ Age, data=survey) EZR pairs() 5 (1) (2) mean µ X µ = N 11

12 X N X = X 1 + X 2 + X X N X X X = n n * 14 median X = n 1( X 1 ) + n 2 ( X 2 ) n n ( X n ) n 1 + n n n (sorting) central tendency R median() Mode R table(x)[which.max(table(x))] (1) (2) (3) (1) (2) * 15 3 (geometric mean) (harmonic mean) Variability 4 (Inter-Quartile Range; IQR) 1/4, 2/4, 3/4 (quartile) 1/4 3/4 25% 75% 2/4 R fivenum() *14 X X X X X C X *15 12

13 Q1, Q2, Q3 50% (Semi Inter-Quartile Range; SIQR) IQR SIQR (variance) V (X µ) 2 V = N * 16 n n 1 (unbiased variance) V ub (X X) 2 V ub = n 1 R var() (standard deviation) R sd() * 17 Mean±2SD * 18 95% EZR EZR wbc OK 2 t EZR wbc OK NA No outliers were identified. 6 *16 *17 *18 2SD

14 Rothman Greenland p * 19 p p a * 20 Welch R t.test(x,y) * 21 b Wilcoxon R wilcox.test(x,y) 2. R prop.test() 6.1 F X Y n X n Y X Y SX<-var(X) SY<-var(Y) SX>SY F0<-SX/SY 1 DFX<-length(X)-1 2 DFY<-length(Y)-1 F 1-pf(F0,DFX,DFY) F0 var.test(x,y) X C C X var.test(x~c) EZR 2 F X 1 C EZR, OK survey Height Sex OK Output *19 fmsb pvalueplot() *20 shapiro.test() Shapiro-Wilk *21 F F 2 t Welch F Welch 14

15 6.2 Welch t t 0 = E(X) E(Y) / S X /n X + S Y /n Y ϕ t ϕ ϕ = (S X /n X + S Y /n Y ) 2 {(S X /n X ) 2 /(n X 1) + (S Y /n Y ) 2 /(n Y 1)} R t.test(x,y,var.equal=f) var.equal Welch t.test(x,y) X C t.test(x~c) survey t.test(height ~ Sex, data=survey) EZR 2 t Height Sex No (Welch test) OK Output p * 22 stripchart() V <- rnorm(100,10,2) W <- rnorm(60,12,3) X <- c(v, W) C <- as.factor(c(rep("v", length(v)), rep("w", length(w)))) x <- data.frame(x, C) x <- stack(list(v=v, W=W)) names(x) <- c("x", "C") * 23 * 24 stripchart(x~c, data=x, method="jitter", vert=true) Mx <- tapply(x$x, x$c, mean) Sx <- tapply(x$x, x$c, sd) Ix <- c(1.1, 2.1) points(ix, Mx, pch=18, cex=2) arrows(ix, Mx-Sx, Ix, Mx+Sx, angle=90, code=3) *22 R barplot() arrows() *23 EZR *24 EZR X C jitter EZR 15

16 6.3 paired-t 0 R X Y paired-t t.test(x,y,paired=t) t.test(x-y,mu=0) survey R t.test(survey$wr.hnd, survey$nw.hnd, paired=true) 1 cm 1 cm 1 cm Diff.Hnd <- survey$wr.hnd - survey$nw.hnd C.Hnd <- ifelse(abs(diff.hnd)<1, 1, ifelse(diff.hnd>0, 2, 3)) matplot(rbind(survey$wr.hnd, survey$nw.hnd), type="l", lty=1, col=c.hnd, xaxt="n") axis(1, 1:2, c("wr.hnd", "NW.Hnd")) EZR 2 paired t Wr.Hnd NW.Hnd [OK] Output 5% 6.4 Wilcoxon Wilcoxon t Mann-Whitney U Kendall S survey (Height) R wilcox.test(height ~ Sex, data=survey) EZR 2 Mann-Whitney U Height Sex OK Output 6.5 Wilcoxon Wilcoxon t survey Wr.Hnd NW.Hnd 5% R 16

17 wilcox.test(survey$wr.hnd, survey$nw.hnd, paired=true) EZR 2 Wilcoxon 1 Wr.Hnd 2 NW.Hnd [OK] n 2 n C % % (Kruskal-Wallis) 5% % 7.1 X, Y, Z * 25 ID (VG) (cm)(height) 1 X X Z Y HEIGHT VG R > sp <- read.delim(" > summary(aov(height ~ VG, data=sp)) *25 R read.delim() 17

18 Df Sum Sq Mean Sq F value Pr(>F) VG ** Residuals Signif. codes: 0 *** ** 0.01 * * Sum Sq VG Sum Sq VG Residuals Sum Sq Mean Sq (Df) VG Mean Sq Residuals Mean Sq F value F Pr(>F) VG 5% EZR sample2.dat sp [ ] [ ] [ URL ] [ :] sp [ URL] [ ] [OK] [OK] ANOVA [ ] [ ] [ one-way ANOVA)] HEIGHT VG [OK] p VG (Bartlett) R Y C bartlett.test(y~c) * 26 p % EZR (Bartlett ) HEIGHT VG [OK] 2 2 Welch R oneway.test() oneway.test(height ~ VG, data=sp) Welch *26 dat bartlett.test(y~c, data=dat) 18

19 > oneway.test(height ~ VG, data=sp) One-way analysis of means (not assuming equal variances) data: HEIGHT and VG F = , num df = 2.00, denom df = 18.77, p-value = EZR aov oneway.test 7.2 (Kruskal-Wallis) Fligner-Killeen (Kruskal-Wallis) R Y C kruskal.test(y~c) Kruskal-Wallis R i (i = 1, 2,..., k; k ) n i N B i B i = n i {R i /n i (N + 1)/2} 2 B = k i=1 B H = 12 B/{N(N + 1)} H A H H H = B i H 1 A N(N 2 1) H H k 4 4 k = 3 5 H H k 1 R kruskal.test(height ~ VG, data=sp) EZR 3 (Kruskal-Wallis ) VG HEIGHT [OK] Fligner-Killeen Bartlett R fligner.test(height ~ VG, data=sp) EZR 19

20 7.3 (Bonferroni) (Holm) (Scheffé) (Tukey) HSD (Dunnett) (Williams) FDR(False Discovery Rate) HSD * 27 HSD FDR R pairwise.t.test() pairwise.wilcox.test() pairwise.prop.test() FDR fmsb pairwise.fisher.test() Fisher Bonferroni p p.adjust.method="fdr" p p p p Benjamini Y, Hochberg Y: Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. Royal Stat. Soc. B, 57: , Bonferroni Holm FDR k p p(1)<p(2)<...<p(k) Bonferroni p(1) /k p(i) /k Holm p(i) /i fdr p(k) p(k-1) (k-1)/k p p(i)< i/k i R pairwise.*.test() Bonferroni p k Holm i p i fdr i p k/i p FDR Rcmdr EZR R pairwise.t.test(sp$height, sp$vg, p.adjust.method="bonferroni") 2 p * 28 pairwise.wilcox.test(sp$height,sp$vg,p.adjust.method="bonferroni") p.adjust.method p.adjust.method="holm" FDR p.adjust.method="fdr" R R accept *27 t *28 "bonferroni" "bon" pairwise.* data= attach() 20

21 TukeyHSD(aov(HEIGHT ~ VG, data=sp)) HSD EZR 3 one-way ANOVA) HEIGHT VG 2 (post-hoc ) OK TukeyHSD 2 95% Tukey p 5% Z Y > TukeyHSD(AnovaModel.3, "factor(vg)") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = HEIGHT ~ factor(vg), data = sp, na.action = na.omit) $ factor(vg) diff lwr upr p adj Y-X Z-X Z-Y Dunnett Dunnett (mmhg ) 5, 8, 3, 10, , 12, 30, 16, , 25, 17, 40, 23 Dunnett R bpdown Dunnett bpdown <- data.frame( medicine=factor(c(rep(1,5),rep(2,5),rep(3,5)), labels=c(" "," "," ")), sbpchange=c(5, 8, 3, 10, 15, 20, 12, 30, 16, 24, 31, 25, 17, 40, 23)) summary(res1 <- aov(sbpchange ~ medicine, data=bpdown)) library(multcomp) res2 <- glht(res1, linfct = mcp(medicine = "Dunnett")) confint(res2, level=0.95) summary(res2) multcomp glht() linfct Dunnett multcomp simtest()

22 EZR URL bpdown URL OK URL OK 3 one-way ANOVA) sbpchange medicine 2 (Dunnett ) OK medicine o.placebo 1.usual 2.newdrug Dunnett 3 (Kruskal- Wallis ) sbpchange medicine 2 (post-hoc Steel ) OK Steel 8 MASS survey R require(mass) MASS plot(wr.hnd ~ Height, data=survey) pch=as.integer(sex) EZR MASS survey OK x Height y Wr.Hnd 2 [OK] Sex [OK] [OK] (correlation) (positive correlation) (negative correlation) 22

23 (apparent correlation) * 29 (spurious correlation) * 30 (Pearson s Product Moment Correlation Coefficient) r X Y X Y [ 1, 1] r = 1 r = 1 2 r = 0 X X Y Ȳ r = ni=1 (X i X)(Y i Ȳ) n i=1 (X i X) 2 n i=1 (Y i Ȳ) 2 p p 5% p = % p t 0 = r n 2 1 r 2 n 2 t survey R methos=spearman cor.test(survey$height, survey$wr.hnd) EZR Pearson Height Wr.Hnd Ctrl OK Rcmdr Pearson s product-moment correlation data: survey$height and survey$wr.hnd t = , df = 206, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: sample estimates: cor *29 *

24 r = % [0.50, 0.69] * 31 p-value < 2.2e ρ a X i R i Y i Q i ρ = 1 6 n(n 2 1) n (R i Q i ) 2 10 i=1 T = ρ n 2 1 ρ 2 n 2 t τ τ = (A B) n(n 1)/2 A B R cor.test() method="spearman" method="kendall" EZR Spearman Spearman Kendall OK a ρ r s % (A) (B) (0, 1, 2, 5, 10 µg/l) (0.24, 0.33, 0.54, 0.83, 1.32) y x y = bx + a a b a b f (a, b) = 5 (y i bx i a) 2 i=1 f (a, b) a b *31 95% 24

25 b = x i y i /5 x i /5 y i /5 i=1 i=1 i= x 2 i /5 x i /5 i=1 i=1 a = 5 5 y i /5 b x i /5 i=1 i=1 a b 0.67 * 32 R lm() linear model y <- c(0.24, 0.33, 0.54, 0.83, 1.32) x <- c(0, 1, 2, 5, 10) # res <- lm(y ~ x) # summary(res) # plot(y ~ x) abline(res) # 0.67 ( res$coef[1])/res$coef[2] Call: lm(formula = y ~ x) Residuals: Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) ** x *** --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 3 degrees of freedom Multiple R-squared: , Adjusted R-squared: F-statistic: 316 on 1 and 3 DF, p-value: a = b = % *32 25

26 (0.9875) Adjusted R-squared p-value µg/l EZR workingcurve [OK] [var1] y numeric Enter [var2] [x] x x y y [OK] y x [OK] VIF survey * 33 R res <- lm(wr.hnd ~ Height, data=survey) summary(res) EZR : survey survey Wr.Hnd Height [OK] *33 survey MASS 26

27 8.4 b a t Y X Y = a 0 + b 0 X + e e 0 σ 2 a a 0 (σ 2 /n)(1 + M 2 /V) M V x Q σ 2 Q/σ 2 (n 2) t 0 (a 0 ) = n(n 2)(a a0 ) (1 + M2 /V)Q (n 2) t a 0 a 0 0 a 0 = 0 t 0 (0) t t 0 (a 0 ) (n 2) t 95% t 97.5% 2 t 0 t 0 t 0 (b) = n(n 2)Vb Q (n 2) t 0 R EZR Pr(> t ) 9 y (saturate) airquality Ozone ppb Solar.R 8:00 12: Langley Wind LaGuardia 7:00 10:00 Temp Month Day R plot(ozone ~ Solar.R, data=airquality) res <- lm(ozone ~ Solar.R, data=airquality) abline(res) 27

28 summary(res) EZR datasets airquality OK airquality x Solar.R y Ozone [OK] Ozone Solar.R OK R EZR Call: lm(formula = Ozone ~ Solar.R, data = airquality) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) ** Solar.R *** --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 109 degrees of freedom (42 observations deleted due to missingness) Multiple R-Squared: ,Adjusted R-squared: F-statistic: on 1 and 109 DF, p-value: Ozone = Solar.R F p Adjusted R-squared % 10 (2011) Lesson ( (2009) R. [ ] [ ] (1983).. (2003) R.. (2007) R.. (2012) EZR EBM * 34. *