R EZR 2013 11 5 *1 1 R 2 1.1 R [2013 11 5 ]................................ 2 1.2 R................................................ 3 1.3 Rgui......................................... 3 1.4 EZR................................................... 4 2 5 2.1 -.......................................... 5 2.2..................................... 6 2.3 (cross-over design)..................................... 6 3 8 4 9 4.1.............................................. 10 4.2................................................. 10 5 11 6 13 6.1 F.......................................... 14 6.2 Welch t........................................... 15 6.3..................................... 16 6.4 Wilcoxon............................................. 16 6.5 Wilcoxon......................................... 16 7 3 17 7.1................................................ 17 7.2 (Kruskal-Wallis) Fligner-Killeen.................. 19 7.3....................................... 20 7.4 Dunnett.............................................. 21 *1 http://minato.sip21c.org/ebhc-text.pdf 1

8 22 8.1................................................ 22 8.2..................................................... 22 8.3............................................. 24 8.4...................................... 27 9 27 10 28 e-mail: minato-nakazawa@umin.net 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 [2013 11 5 ] Windows CRAN R-3.0.2 *6 *7 SDI (separate windows) Macintosh Mac OS X CRAN R-3.0.2.pkg Tcl/Tk tools *8 *2 R http://www.okada.jp.org/rwiki/ *3 http://cran.md.tsukuba.ac.jp/ *4 http://essrc.hyogo-u.ac.jp/cran/ *5 http://cran.ism.ac.jp/ *6 R-3.0.2-win.exe *7 R 1 MDI Rcmdr/EZR *8 http://aoki2.si.gunma-u.ac.jp/r/begin.html 2

Linux Debian RedHat ubuntu CRAN CRAN R-3.0.2.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() <- 1 4 6 3 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

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 http://www.jichi.ac.jp/saitama-sct/saitamahp.files/statmed.html *12 2013 8 23 Rcmdr 2.0.0 RcmdrPlugin.EZR EZR R Windows R EZR 4

: 3 3 triplicate 3 3 2 2 3 2 2 duplicate triplicate triplicate 96 3 3 2 2 60 40 3 2 2 0.1% 5% 3 3 2 R.A. Fisher - 2.1 - t 2 t t (RA) 5

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): 902-906. 574 43 531 2 1 266 2 265 1 A 34 2 B 15 2 2 250 A 1 232 B 45 29 2 408 * 13 PID Ord.LD LD DL LD L D T0.Alb T1.Alb *13 http://minato.sip21c.org/pubhealthpractice/crossovertest.txt web 6

PID Ord.LD LD T0.Alb T1.Alb 1 LD L 15 21 2 LD L 17 19 3 LD L 19 25 4 LD L 21 23 5 LD L 23 27 6 DL D 15 34 7 DL D 17 27 8 DL D 19 39 9 DL D 21 28 10 DL D 23 41 1 LD D 15 35 2 LD D 17 28 3 LD D 19 40 4 LD D 21 29 5 LD D 23 42 6 DL L 15 22 7 DL L 17 20 8 DL L 19 26 9 DL L 21 24 10 DL L 23 28 2 2 ) 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) 22944.1 1 402.8 16 911.3843 1.554e-15 *** Factor1.LD 291.6 1 402.8 16 11.5829 0.003634 ** Factor2.Ord.LD 0.0 1 402.8 16 0.0000 1.000000 Factor1.LD:Factor2.Ord.LD 2.5 1 402.8 16 0.0993 0.756738 Time 980.1 1 154.8 16 101.3023 2.510e-08 *** Factor1.LD:Time 291.6 1 154.8 16 30.1395 4.942e-05 *** Factor2.Ord.LD:Time 0.0 1 154.8 16 0.0000 1.000000 Factor1.LD:Factor2.Ord.LD:Time 2.5 1 154.8 16 0.2584 0.618160 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 Pr(>F) Factor1.LD 5% Time 5% Factor1.LD:Time 5% Factor2.Ord.LD Pr(>F) 2 7

3 60 kg, 66 kg, 75 kg R mean(c(60,66,75)) (60+66+75)/3 Microsoft Excel 10 ID (cm) (kg) 1 170 70 2 172 80 3 166 72 4 170 75 5 174 55 6 199 92 7 168 80 8 183 78 9 177 87 10 185 100 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

Dataset <- read.delim("desample.txt") Rcmdr a Dataset OK Excel OK Windows R-2.9.0 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

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) 10 20 10 hist(survey$age, breaks=1:8*10, right=false) Rcmdr EZR survey Age 70 16.75 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

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 3 +1.5x 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

X N X = X 1 + X 2 + X 3 +... + 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 2 +... + 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

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 97.5 1.959964... 2 2 13

Rothman Greenland p * 19 p p 0.05 0.01 1. 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

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 2013 8 EZR 15

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

wilcox.test(survey$wr.hnd, survey$nw.hnd, paired=true) EZR 2 Wilcoxon 1 Wr.Hnd 2 NW.Hnd [OK] 7 3 3 2 n 2 n C 2 1 1 5% 3 1 1 5% (Kruskal-Wallis) 5% 2 1 1 5% 7.1 X, Y, Z * 25 ID (VG) (cm)(height) 1 X 161.5 2 X 167.0 22 Z 166.0 37 Y 155.5 HEIGHT VG R > sp <- read.delim("http://minato.sip21c.org/grad/sample2.dat") > summary(aov(height ~ VG, data=sp)) *25 http://minato.sip21c.org/grad/sample2.dat R read.delim() 17

Df Sum Sq Mean Sq F value Pr(>F) VG 2 422.72 211.36 5.7777 0.006918 ** Residuals 34 1243.80 36.58 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 * Sum Sq VG Sum Sq 422.72 VG Residuals Sum Sq 1243.80 Mean Sq (Df) VG Mean Sq 211.36 Residuals Mean Sq 36.58 F value 1 2 2 34 F Pr(>F) 0.006918 VG 5% EZR sample2.dat sp [ ] [ ] [ URL ] [ :] sp [ URL] [ ] [OK] http://minato.sip21c.org/grad/sample2.dat [OK] ANOVA [ ] [ ] [ one-way ANOVA)] HEIGHT VG [OK] p VG (Bartlett) R Y C bartlett.test(y~c) * 26 p 0.5785 5% 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

> oneway.test(height ~ VG, data=sp) One-way analysis of means (not assuming equal variances) data: HEIGHT and VG F = 7.5163, num df = 2.00, denom df = 18.77, p-value = 0.004002 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

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: 289-300, 1995. 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

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 -2.538889-8.3843982 3.30662 0.5423397 Z-X 5.850000-0.9598123 12.65981 0.1038094 Z-Y 8.388889 2.3382119 14.43957 0.0048525 7.4 Dunnett Dunnett 5 3 1 (mmhg ) 5, 8, 3, 10, 15 1 20, 12, 30, 16, 24 1 31, 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 0.993 simtest() 0.994 21

EZR URL bpdown URL OK URL http://minato.sip21c.org/bpdown.txt 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] 8.1 8.2 2 (correlation) (positive correlation) (negative correlation) 22

(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 = 0.034 5% 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 = 10.7923, df = 206, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.5063486 0.6813271 sample estimates: cor 0.6009909 *29 *30 15 23

r = 0.60 95% [0.50, 0.69] * 31 p-value < 2.2e-16 2.2 10 16 0.7 0.4 0.7 0.2 0.4 ρ 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 8.3 98% (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

b = 5 5 5 x i y i /5 x i /5 y i /5 i=1 i=1 i=1 2 5 5 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 (0.67 - res$coef[1])/res$coef[2] Call: lm(formula = y ~ x) Residuals: 1 2 3 4 5-0.02417-0.04190 0.06037 0.02718-0.02147 Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) 0.26417 0.03090 8.549 0.003363 ** x 0.10773 0.00606 17.776 0.000388 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 Residual standard error: 0.04894 on 3 degrees of freedom Multiple R-squared: 0.9906, Adjusted R-squared: 0.9875 F-statistic: 316 on 1 and 3 DF, p-value: 0.0003882 a = 0.26417 b = 0.10773 98.75% *32 25

(0.9875) Adjusted R-squared p-value 0.67 3.767084 3.8 µ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

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 1973 5 1 9 30 154 Ozone ppb Solar.R 8:00 12:00 4000 7700 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

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 -48.292-21.361-8.864 16.373 119.136 Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) 18.59873 6.74790 2.756 0.006856 ** Solar.R 0.12717 0.03278 3.880 0.000179 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 Residual standard error: 31.33 on 109 degrees of freedom (42 observations deleted due to missingness) Multiple R-Squared: 0.1213,Adjusted R-squared: 0.1133 F-statistic: 15.05 on 1 and 109 DF, p-value: 0.0001793 Ozone = 18.599 + 0.127 Solar.R F p 0.0001793 Adjusted R-squared 0.11 10% 10 (2011) Lesson 3. 2937 (http://www.igaku-shoin.co.jp/paperdetail.do?id=pa02937_06) (2009) R. [ ] [ ] (1983).. (2003) R.. (2007) R.. (2012) EZR EBM * 34. *34 http://www.jichi.ac.jp/saitama-sct/saitamahp.files/statmed.html 28