win.metafile() emf Microsoft PowerPoint OpenOffice.org Draw 2 R CRAN R (2004) [ ](2004) (2005) The R Tips R The R Tips R 6

Transcription

1 R nminato@med.gunma-u.ac.jp 1 R R GNU General Public License Version 2 1 GIS 2 R Fisher iris Gehan MASS Windows gehan Windows Macintosh Linux FreeBSD R CDC EPIINFO EPIINFO Windows R R SPSS R 3 R S test.r source("test.r") R R Commander Rcmdr pdf postscript png jpeg Windows (emf) emf Windows R 1 LGPL (Lesser GNU General Public License) Version CRAN ( CRAN R built-in base, datasets, grdevices, graphics, grid, methods, splines, stats, stats4, tcltk, tools, utils Recommended built-in KernSmooth, MASS, boot, class, cluster, foreign, lattice, mgcv, nlme, nnet, rpart, spatial, survival Windows library(survival) require(survival) search().packages(all.avail=t) detach(package:survival) Nature R Morris RJ, Lewis OT, Godfray HCJ: Experimental evidence for apparent competition in a tropical forest food web. Nature 428: ,

2 win.metafile() emf Microsoft PowerPoint OpenOffice.org Draw 2 R CRAN R (2004) [ ](2004) (2005) The R Tips R The R Tips R 6 RjpWiki 7 R 2 R Windows Linux Mac OS X \ Rgui R 8.Rprofile.RData > > R + R Source source( ) Windows / \\ R bin Windows 2000/XP R R R 2.1 R 4 ftp://ftp.u-aizu.ac.jp/pub/lang/r/cran/ (S) R_USER.Renviron R_USER="c:/work" proxy Windows proxy --internet2 2

3 > q() R 9 X x " NA comment(x)<-"test" > x <- 7 > x [1] 7 > comment(x) <- " " > x [1] 7 > comment(x) [1] " " > 6 -> x > x [1] 6 > comment(x) NULL > names(x) <- " " > x 6 > x[1] <- 4 > x 4 > z <<- 5 # R html pdf Html t t.test() > help(t.test); #?t.test Fisher fisher.test(stats) Fihser 9 3

4 > help.search("fisher") R > example(barplot) R c() : > x <- c(2,7,11:19) > x [1] R function() list() > x <- 2 > z <- function(a) { x <<- x+a } > print(z(5)) [1] 7 > x [1] 7 > meansd <- function(x) {list(mean=mean(x),sd=sd(x))} > meansd(c(1,5,8)) $mean [1] $sd [1] CRAN CRAN > options(repos=" R.Rprofile install.packages(" ") vcd dep=t dependency TRUE > install.packages("vcd",dep=t) 4

5 2.2 R factor character ordered numeric integer logical data.frame list ts matrix vector data.frame numeric double single C Fortran complex dat$c <- as.factor(dat$c) dat$s <- as.character(dat$s) dat$i <- as.ordered(dat$i) dat$x <- as.numeric(dat$x) str(dat) names(dat) R 155 cm 160 cm 170 cm c(155,160,170) seq() rep() mean(dat) sd(dat) cm 50 kg A 160 cm 55 kg B 170 cm 70kg C data.frame(ht=c(155,160,170),wt=c(50,55,70)) $ attach() attach() detach() > dat <- c(155,160,170) > dat2 <- seq(155,170,by=5) > dat3 <- rep(160,3) > mean(dat) [1] > sd(dat) [1] > dat4 <- data.frame(ht=c(155,160,170),wt=c(50,55,70)) > dat4$ht [1] > attach(dat4) > wt [1] > detach(dat4)

6 dat > dat <- matrix(c(20,12,10,18),nc=2) > rownames(dat) <- c(, ) > colnames(dat) <- c(, ) 2 3 matrix() 1 2 nc=2 2 1 chisq.test(dat) ID (cm) (kg) Microsoft Excel OpenOffice.org 11 calc PID, HT, WT Microsoft Excel (F) (T) (*.txt) xls (S) OK Excel R R sample.txt R dat <- read.delim("sample.txt") dat R de() 11 Microsoft Office 12 Windows Excel calc R dat <- read.delim("clipboard") 6

7 > attach(dat) > cor.test(ht,wt) > detach(dat) 3.3 (NA) R NA Excel 13 Excel R dat 1 datc > datc <- subset(dat,complete.cases(dat),drop=t) (recovery rate) 80% 80/100) (effective recovery rate) 75% 75/100 80% 3.4 Microsoft Access Oracle postgresql PHP4 Apache httpd 14 Tcl/Tk R RODBC R Dbase Oracle postgresql

8 4 R 4.1 ID (cm) (kg) M F M M F F F M F M 45 1 PID, HT, WT, SEX, AGE R sample.txt > dat <- read.delim("sample.txt") dat c() matrix nc=5 5 1 byrow=t as.data.frame() colnames() 1 Factor M read.delim() > dat <- as.data.frame(matrix(c( + 1, 170, 70, 1, 54, + 2, 162, 50, 2, 34, + 3, 166, 72, 1, 62, + 4, 170, 75, 1, 41, + 5, 164, 55, 2, 37, + 6, 159, 62, 2, 55, + 7, 168, 80, 2, 67, + 8, 183, 78, 1, 47, + 9, 157, 47, 2, 49, + 10, 185, 100, 1, 45),nc=5,byrow=T)) + colnames(dat) <- c( PID, HT, WT, SEX, AGE ) > dat$sex <- as.factor(dat$sex) > levels(dat$sex) <- c( M, F ) data.frame() 8

9 4.2 R [ ] > attach(dat) > mean(ht[sex== M ]) [1] > detach(dat) [dat$sex== M ] > cnms <- function(x,c) { list(n=nrow(x[c]),mean=mean(x[c]),sd=sd(x[c])) } N mean sd cnms() > attach(dat) > cnms(ht,sex== M ) $N [1] 5 $mean [1] $sd [1] > detach(dat) == is.na() 40! > attach(dat) > overforty <- AGE>=40 > cnms(ht,overforty) $N [1] 8 $mean [1] $sd [1] > detach(dat) & 40 9

10 > attach(dat) > males <- SEX== M > cnms(wt,overforty&males) $N [1] 5 $mean [1] 79 $sd [1] > detach(dat) [ ] 4.3 R tapply() > attach(dat) > tapply(ht,sex,mean) M F > detach(dat) 5 TFR PowerPoint OpenOffice.org Impress PowerPoint 16 OpenOffice.org Draw Adobe Illustrator 10

11 (2005) (ISBN ) X R barplot(table(x)) M ± table(x) c() names() barplot() OpenOffice.org Draw 17 table() 11

12 > ob <- c(4,1,2,12,97) > names(ob) <- c("+++","++","+"," ","-") > barplot(ob,ylim=c(0,100),main=" \n ") 2 barplot() > ob <- c(4,1,2,12,97) > names(ob) <- c("+++","++","+"," ","-") > ii <- barplot(matrix(ob,nrow(ob)),beside=f,ylim=c(0,120), + main=" ") > oc <- ob > for (i in 1:length(ob)) { oc[i] <- sum(ob[1:i])-ob[i]/2 } > text(ii,oc,paste(names(ob))) ± ± 7 50 R > obm <- c(0,0,1,5,47) > obf <- c(4,1,1,7,50) > obx <- cbind(obm,obf) > rownames(obx) <- c("+++","++","+"," ","-") > colnames(obx) <- c(" "," ") > ii <- barplot(obx,beside=f,ylim=c(0,70),main=" ") > oc <- obx > for (i in 1:length(obx[,1])) { oc[i,1] <- sum(obx[1:i,1])-obx[i,1]/2 } > for (i in 1:length(obx[,2])) { oc[i,2] <- sum(obx[1:i,2])-obx[i,2]/2 } > text(ii[1],oc[,1],paste(rownames(obx))) > text(ii[2],oc[,2],paste(rownames(obx))) 12

13 100% R 2 (%) 4 horiz=t > ob <- c(4,1,2,12,97) > obp <- ob/sum(ob)*100 > names(obp) <- c("+++","++","+"," ","-") > ii <- barplot(matrix(obp,nrow(obp)),horiz=t,beside=f,xlim=c(0,100), + xlab="(%)",main=" ") > oc <- obp > for (i in 1:length(obp)) { oc[i] <- sum(obp[1:i])-obp[i]/2 } > text(oc,ii,paste(names(obp))) barplot() dotchart() 13

14 > obm <- c(0,0,1,5,47) > obf <- c(4,1,1,7,50) > obx <- cbind(obm,obf) > rownames(obx) <- c("+++","++","+"," ","-") > colnames(obx) <- c(" "," ") > dotchart(obx) > dotchart(t(obx)) 100% 18 R pie() 19 R > ob <- c(4,1,2,12,97) > names(ob) <- c("+++","++","+"," ","-") > pie(ob) 5.2 Excel R hist() Cleveland Cleveland WS (1985) The elements of graphing data. Wadsworth, Monterey, CA, USA. p.264 R help 19 R-1.5 piechart()

15 > attach(dat) > hist(ht,main=" ") > detach(dat) qqnorm() qqline() > qqnorm(dat$ht,main=" ",ylab=" (cm)") > qqline(dat$ht,lty=2) stem and leaf plot 5 10 R stem() > stem(dat$ht) box and whisker plot boxplot (median) 1/4 (first quartile) 1/4 (third quartile) 1.5 R boxplot() > boxplot(dat$ht) R stripchart() vert=t > attach(dat) > mht <- tapply(ht,sex,mean) > sht <- tapply(ht,sex,sd) > IS <- c(1,2)+0.15 > stripchart(ht~sex,method="jitter",vert=t,ylab=" (cm)") > points(is,mht,pch=18) > arrows(is,mht-sht,is,mht+sht,code=3,angle=90,length=.1) > detach(dat) 21 source(" stem() gstem() 15

16 R plot() plot() pch points() symbols() 22 matplot() matpoints() pairs() text() identify() > plot(dat$ht,dat$wt,pch=paste(dat$sex),xlab=" (cm)",ylab=" (kg)") R stars() 5.3 maptools ESRI GIS 23 6 N(µ, σ 2 ) µ σ 2 µ σ 2 central tendency variability parameter (location parameter) (scale

17 6.1 1 table(c) Q1 Q3 fivenum(x) NROW(X) length(x) sum(x) max(x) min(x) mean(x) sum(x)/length(x) exp(mean(log(x))) prod(x)^(1/nrow(x)) 1/(mean(1/X)) median(x) fivenum(x)[4]-fivenum(x)[2] (fivenum(x)[4]-fivenum(x)[2])/2 var(x) sum((x-mean(x))^2)/(length(x)-1) sd(x) sqrt(var(x)) table(c1,c2) cor(x,y) cor(x,y,method="spearman") 6.3 R try(data()) ChickWeight Chick Diet Time weight X > data(chickweight) > attach(chickweight) > X <- weight[time==20] > detach(chickweight) R d p q r 15 t 97.5% qt(0.975,df=15) rnorm(100,0,1) 17

18 norm t t F f chisq wilcox help.search() curve() 10 2 (0,20) > curve(dnorm(x,10,2),0,20) R p.value 7.3 Shapiro-Wilk Shapiro-Wilk Z i = (X i µ)/σ Z i X N(0, 1) c(i) = E[Z(i)], d ij = Cov(Z(i), Z(j)) X(1) < X(2) <... < X(n) c(1), c(2),...c(n) σ n ˆ(σ) = i=1 a ix(i) S 2 = n i=1 (X i X) 2 W = (kˆσ 2 )/S 2 k n i=1 (ka i) 2 = 1 X > shapiro.test(x) X R length(x) W 24 5% 1% 18

19 7.4 1 shapiro.test(x) t.test(x,mu= ) binom.test(table(b)[2],length(b),p= ) 7.5 X Y µ X µ Y µ X = mean(x) = X/n µ Y = mean(y ) = Y/n H0 : µ X = µ Y H1 : µ X µ Y H1 µ X > µ Y µ X < µ Y t 0 t 0 5% t 0 t 2.5% t 0 t 2.5% 97.5% t t 0 t 25 X Y X Y H0 : µ X µ Y H1 : µ X < µ Y t 0 5% t 0 t 5% 95% R t.test() alt="greater" alt="less" X Y n X n Y 26 V z 0 = E(X) E(Y ) / V/n X + V/n Y F 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) 2. Welch 25 t t 0 t t 0 1 R 1-pt(t0, ) [11] 27 Mann-Whitney U Wilcoxon 19

20 S S<-(DFX*SX+DFY*SY)/(DFX+DFY) t0<-abs(mean(x)-mean(y))/sqrt(s/length(x)+s/length(y)) DFX+DFY t X Y (1-pt(t0,DFX+DFY))*2 R t.test(x,y,var.equal=t) F t.test(x~c,var.equal=t) t.test(x,y,var.equal=t,alternative="less") alternative="less" X<Y X>=Y Welch 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) t.test(x~c) Fisher setosa virginica Sepal.Length t Fisher versicolor 3 3 setosa virginica p-value > data(iris) > setosa.sl <- iris$sepal.length[iris$species=="setosa"] > virginica.sl <- iris$sepal.length[iris$species=="virginica"] > vareq <- var.test(setosa.sl,virginica.sl)$p.value >= 0.05 > t.test(setosa.sl,virginica.sl,var.equal=vareq) Welch Two Sample t-test data: setosa.sl and virginica.sl t = , df = , p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: sample estimates: mean of x mean of y setosa virginica virginica 50 Sepal.Length Petal.Length 2 paired-t 0 R X Y paired-t t.test(x,y,paired=t) t.test(x-y,mu=0) 20

21 R t p-value virginica > data(iris) > virginica <- subset(iris,species=="virginica",drop=t) > attach(virginica) > t.test(sepal.length,petal.length,paired=t) Paired t-test data: Sepal.Length and Petal.Length t = , df = 49, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: sample estimates: mean of the differences > wilcox.test(sepal.length,petal.length,paired=t) Wilcoxon signed rank test with continuity correction data: Sepal.Length and Petal.Length V = 1275, p-value = 7.481e-10 alternative hypothesis: true mu is not equal to 0 > detach(virginica) n 2 nc % % (Kruskal-Wallis) 5% 2 1 5% 2 (1997),

22 post hoc 8.2 R chickwts 71 R Console?chickwts Anonymous (1948) Biometrika, 35: (g) (casein) 368, 390, 379, 260, 404, 318, 352, 359, 216, 222, 283, 332 (horsebean) 179, 160, 136, 227, 217, 168, 108, 124, 143, 140 (linseed) 309, 229, 181, 141, 260, 203, 148, 169, 213, 257, 244, 271 (meatmeal) 325, 257, 303, 315, 380, 153, 263, 242, 206, 344, 258 (soybean) 243, 230, 248, 327, 329, 250, 193, 271, 316, 267, 199, 171, 158, 248 (sunflower) 423, 340, 392, 339, 341, 226, 320, 295, 334, 322, 297, 318 chickwts weight feed > data(chickwts) > attach(chickwts) > layout(cbind(1,2)) > boxplot(weight~feed,ylab=" (g)") > stripchart(weight~feed,vert=t,method="jitter",ylab=" (g)") > detach(chickwts) weight feed

23 R Console summary(aov(weight~feed)) anova(lm(weight~feed)) > attach(chickwts) > summary(aov(weight~feed)) Df Sum Sq Mean Sq F value Pr(>F) feed e-10 *** Residuals Signif. codes: 0 *** ** 0.01 * > detach(chickwts) * Pr(>F) Sum Sq (sum of squares) feed Sum Sq feed Residuals Sum Sq Mean Sq (mean square) (Df) feed Mean Sq Residuals Mean Sq 3009 F value F F > 1 1-pf(15.365,5,65) Pr(>F) 5.936e feed 5% 6 (Bartlett) R bartlett.test( ~ ) kruskal.test() bartlett.test(weight~feed) p-value % 29 1-pf(15.365,5,65) 3 1 F value 23

24 > attach(chickwts) > print(res.bt <- bartlett.test(weight~feed)) Bartlett test of homogeneity of variances data: weight by feed Bartlett s K-squared = , df = 5, p-value = 0.66 > ifelse(res.bt$p.value<0.05, + cat(" Bartlett p=",res.bt$p.value,"\n"), + summary(aov(weight~feed))) [[1]] Df Sum Sq Mean Sq F value Pr(>F) feed e-10 *** Residuals Signif. codes: 0 *** ** 0.01 * > detach(chickwts) 8.3 R R Holm R > attach(chickwts) > pairwise.t.test(weight,feed) Pairwise comparisons using t tests with pooled SD data: weight and feed casein horsebean linseed meatmeal soybean horsebean 2.9e linseed meatmeal e soybean sunflower e e P value adjustment method: holm > detach(chickwts) Fisher 3 (Petal.Width) 24

25 (case control study) R 2 table() mosaicplot() > pid <- 1:13 > sex <- as.factor(c(rep(1,6),rep(2,7))) > levels(sex) <- c(" "," ") > disease <- as.factor(c(1,1,1,2,2,2,1,1,1,1,2,2,2)) > levels(disease) <- c(" "," ") > print(ctab <- table(sex,disease)) disease sex > mosaicplot(ctab,main="2 2 ") chisq.test(ctab) (cohort study) 25

26 a a /200 = /200 = /200 = /200 = 32.5 χ 2 c = (80 68) 2 / (55 67) 2 / (20 32) 2 / (45 33) 2 /32.5 = pchisq(13.128,1) % R > X <- matrix(c(80,20,55,45),nr=2) > chisq.test(x) Pearson s Chi-squared test with Yates continuity correction data: X X-squared = , df = 1, p-value = > smoker <- c(80,55) > pop <- c(100,100) > prop.test(smoker,pop) Fisher X fisher.test(x) % 2 2 fisher.test() 95%

27 (1) (prevalence) point prevalence prevalence prevalence (cumulative incidence) (risk) 20 (incidence rate) International Epidemiological Association Last JM [Ed.] A Dictionary of Epidemiology, 4th Ed. (Oxford Univ. Press, 2001) incidence (odds) (disease-odds) (exposure-odds) (2) Relative Risk (risk ratio) cumulative incidence rate ratio (incidence rate ratio) (mortality rate ratio) rate ratio (odds ratio) 2 (attributable risk) (risk difference) (incidence rate difference) (excess risk) (attributable proportion) 1 95% 1 5% 32 (prospective study) (cohort 32 Rothman Greenland 95% Rothman p-value 27

28 (follow-up study) case-control study 33 cross-sectional study 34 (exposure-odds ratio) (disease-odds ratio) a b m 1 c d m 2 n 1 n 2 N a/m 1 c/m 2 = am 2 cm 1 a/b c/d = ad bc a/c b/d = ad bc R fisher.test() ad/bc fisher.test() vcd oddsratio() log=f (a+0.5)(d+0.5) 0 (b+0.5)(c+0.5) summary(oddsratio()) confint(oddsratio())

29 (4/100000)/(2/100000) = 2 ( )/( ) Yule Q test-retest-reliability κ 95% m 1 m 2 X Y X m 1 X m 1 Y m 2 Y m 2 X + Y N X Y N π 1 = X/m 1 π 2 = Y/m 2 RR = π 1 /π 2 = (Xm 2 )/(Y m 1 ) N Bailey 95% RR exp( qnorm(0.975) 1/X 1/m 1 + 1/Y 1/m 2 ) (1) 35 95% 29

30 RR exp(qnorm(0.975) 1/X 1/m 1 + 1/Y 1/m 2 ) (2) RR 2 95% (0.37, 10.9) 2 95% 1 5% > riskratio2 <- function(x,y,m1,m2) { + data <- matrix(c(x,y,m1-x,m2-y,m1,m2),nr=2) + colnames(data) <- c(" "," ") + rownames(data) <- c(" "," ") + print(data) + RR <- (X/m1)/(Y/m2) + RRL <- RR*exp(-qnorm(0.975)*sqrt(1/X-1/m1+1/Y-1/m2)) + RRU <- RR*exp(qnorm(0.975)*sqrt(1/X-1/m1+1/Y-1/m2)) + cat(" : ",RR," 95% = [ ",RRL,", ",RRU," ]\n") + } > riskratio2(4,2,100000,100000) epitools riskratio() rateratio() midp wald boot method="wald" median-unbiased epitools 2 95% (0.37, 10.9) median-unbiased > library(epitools) > rateratio(c(2,4,5*100000,5*100000),method="wald") a, b, c, d OR OR = (ad)/(bc) Cornfield (1956) 95% OR exp( qnorm(0.975) 1/a + 1/b + 1/c + 1/d) OR exp(qnorm(0.975) 1/a + 1/b + 1/c + 1/d) Cornfield Newton % (0.37, 10.9) 30

31 > oddsratio2 <- function(a,b,c,d) { + data <- matrix(c(a,b,a+b,c,d,c+d,a+c,b+d,a+b+c+d),nr=3) + colnames(data) <- c(" "," ") + rownames(data) <- c(" "," "," ") + print(data) + OR <- (a*d)/(b*c) + ORL <- OR*exp(-qnorm(0.975)*sqrt(1/a+1/b+1/c+1/d)) + ORU <- OR*exp(qnorm(0.975)*sqrt(1/a+1/b+1/c+1/d)) + cat(" : ",OR," 95% = [ ",ORL,", ",ORU," ]\n") + } > oddsratio2(4,2,99996,99998) Windows Fisher workspace vcd oddsratio() > fisher.test(matrix(c(4,2,99996,99998), nr=2),workspace= ) Fisher s Exact Test for Count Data data: matrix(c(4, 2, 99996, 99998), nr = 2) p-value = alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: sample estimates: odds ratio > require(vcd) > oddsratio(matrix(c(4,2,99996,99998), nr=2),log=f) [1] vcd oddsratio() 95% confint() % [0.43, 9.39] > require(vcd) > OR <- oddsratio(matrix(c(4,2,99996,99998),nr=2),log=f) > ORCI <- confint(or) > cat(" : ",OR," 95% = [ ",ORCI[1],", ",ORCI[2]," ]\n") 9.3 π 1 π 2 (π 1 π 2 )/π 1 (π 1 π 2 )/(1 π 2 ) π = (X + Y )/(m 1 + m 2 ) (π π 2 )/π Q 1 1 Q = (OR 1)/(OR + 1) 0 31

32 (φ) 1,0 θ 1, θ 2 φ = (π 1 π 2 )(θ 1 θ 2 ) k m χ 2 0 n χ 2 0 /n k m t 0 t 1 C C = φ 2 /(1 + φ 2 ) 0 (t 1)/t V V = φ/ t V vcd assocstats() > require(vcd) > assocstats(matrix(c(4,2,99996,99998),nr=2)) 9.4 κ test-retest reliability vcd agreementplot() a b m 1 c d m 2 n 1 n 2 N P e = (n 1 m 1 /N + n 2 m 2 /N)/N P o = (a + d)/n κ = (P o P e )/(1 P e ) κ κ V (κ) V (κ) = P e /(N (1 P e )) κ/ V (κ) κ = 0 κ 95% κ 95% R

33 > kappa.test <- function(x) { + x <- as.matrix(x) + a <- x[1,1]; b <- x[1,2]; c <- x[2,1]; d <- x[2,2] + m1 <- a+b; m2 <- c+d; n1 <- a+c; n2 <- b+d; N <- sum(x) + Pe <- (n1*m1/n+n2*m2/n)/n + Po <- (a+d)/n + kappa <- (Po-Pe)/(1-Pe) + sek0 <- sqrt(pe/(n*(1-pe))) + sek <- sqrt(po*(1-po)/(n*(1-pe)^2)) + p.value <- 1-pnorm(kappa/seK0) + kappal<-kappa-qnorm(0.975)*sek + kappau<-kappa+qnorm(0.975)*sek + list(kappa=kappa,conf.int=c(kappal,kappau),p.value=p.value) + } > kappa.test(matrix(c(10,3,2,19),nr=2)) $kappa [1] $conf.int [1] $p.value [1] e-05 vcd Kappa() m m κ 36 confint() 37 > require(vcd) > print(mykappa <- Kappa(matrix(c(10,3,2,19),nr=2))) > confint(mykappa) C3 2 C1 C2 mantelhaen.test(c1,c2,c3) TMP <- table(c1,c2,c3) 3 TMP mantelhaen.test(tmp) C 2 2 B1 B2 36 Po Pe weights= weights="equal-spacing" weights="fleiss-cohen" nc 1-(abs(outer(1:nc,1:nc,"-"))/(nc-1))^2 weights="equal-spacing" 1-abs(outer(1:nc,1:nc,"-"))/(nc-1) 2 matrix(c(1,0,0,1),nc=2) 37 κ 95 poor slight fair moderate substantial almost perfect 1 perfect Landis and Koch, 1977, Biometrics, 33: vcd κ = 0 33

34 Woolf vcd woolf.test(table(b1,b2,c)) 2 subset() chisq.test(table(c1,c2)) Fisher fisher.test(table(c1,c2)) library(vcd); summary(oddsratio(table(b1,b2),log=f)) library(vcd); assocstats(table(c1,c2)) library(vcd); Kappa(table(C1,C2)) prop.trend.test(table(b,i)[2,],table(i)) Cochran-Armitage var.test(x~b) F t.test(x~b,var.equal=t) t t.test(x~b) Welch wilcox.test(x~b) Wilcoxon t.test(x,y,paired=t) paired-t bartlett.test(x~c) aov(x~c) C TukeyHSD(aov(X~C)) pairwise.t.test(x,c) Tukey Holm library(multcomp) simtest(x~c,type="dunnett") simtest(x~c,type="williams") kruskal.test(x~c) pairwise.wilcox.test(x,c) Wilcoxon cor.test(x,y) cor.test(x,y,method="spearman") cor.test(x,y,method="kendall") 34

35 11 (Generalized Linear Model) Y = β 0 + βx + ε Y 38 X β 0 β ε 39 R glm() lm() CRAN nls() 11.1 R glm() (1) X1 X2 Y Y (2)(1) (3)dat Y Y (4) C1 C2 Y > glm(y ~ X1+X2) > glm(y ~ X1+X2-1) > glm(y ~., data=dat, family="binomial") > glm(y ~ C1+C2+C1:C2) family "gaussian" 2 family (1) lm(y ~ X1+X2) summary(lm()) lm() (4) lm() * (4) C1*C2 (4) anova(lm(y ~ C1*C2)) res <- glm(y ~ X1+X2) plot(residuals(res)) summary(res) AIC(res) AIC step(res) R 35

36 11.2 t (Y ) (X) t Welch 2 Type I Type IV Type I SS Type II Type III Type II Type II SS 2 Type III Type II Type IV SAS SAS MANOVA Type II Type I Type II 16 IML Type III R anova() aov() Type I anova(lm()) aov() car Anova() Type II Type III type= III library(car) help(anova) Anova() Type II SAS Type II Type III car John Fox An R and S-PLUS companion to applied regression. p.140 Type III A B A B A Type III SAS R library(car) Anova(lm(Y~C1*C2)) Type II 36

37 AIC t > res <- lm(y ~ X+Z) > sdd <- c(0,sd(x),sd(z)) > stb <- coef(res)*sdd/sd(y) > print(stb) 11.4 (multicolinearity) Y X1 X2 lm(y ~ X1+X2) X1 X2 X1 X2 X2 (multicolinearity) 2 1 VIF Variance Inflation Factor; VIF 10 Armitage et al DBP SBP centring R MASS lm.ridge() DAAG Maindonald and Braun, 2003 vif() VIF Armitage et al. 37

38 data(airquality) Ozone ppb Solar.R 8:00 12: Langley Wind LaGuardia 7:00 10:00 Temp Month Day LaGuardia Armitage VIF R > data(airquality) > attach(airquality) > res <- lm(ozone ~ Solar.R+Wind+Temp) > VIF <- function(x) { 1/(1-summary(X)$r.squared) } > VIF(lm(Solar.R ~ Wind+Temp)) > VIF(lm(Wind ~ Solar.R+Temp)) > VIF(lm(Temp ~ Solar.R+Wind)) > summary(res) Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) ** Solar.R * Wind e-06 *** Temp e-09 *** --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 107 degrees of freedom Multiple R-Squared: , Adjusted R-squared: F-statistic: on 3 and 107 DF, p-value: < 2.2e-16 > detach(airquality) 3 VIF 10 summary(res) 5% Adjusted R-squared 3 60% 11.5 (1) (2) 38

39 AIC 11.6 res Wind plot(residuals(res)~res$model$wind) plot(residuals(res)) predict() Wind 95% Wind Ozone 42 95% > data(airquality) > attach(airquality) > res <- lm(ozone ~ Solar.R+Wind+Temp) > EW <- seq(min(wind),max(wind),len=100) > ES <- rep(mean(solar.r,na.rm=t),100) > ET <- rep(mean(temp,na.rm=t),100) > Ozone.EWC <- predict(res,list(wind=ew,solar.r=es,temp=et),interval="conf") > plot(ozone~wind) > lines(ew,ozone.ewc[,1],lty=1) > lines(ew,ozone.ewc[,2],lty=2,col="blue") > lines(ew,ozone.ewc[,3],lty=2,col="blue") > detach(airquality) 42 95% 95% predict() interval="conf" interval="pred" 39

40 f(x, θ) {x 1, x 2,..., x n } θ L(θ) = f(x 1, θ)f(x 2, θ) f(x n, θ) L(θ) θ θ L(θ) θ θ θ ln L(θ) µ (1995) -2 λ R loglik() % 2 3 > data(airquality) > attach(airquality) > res.3 <- lm(ozone ~ Solar.R+Wind+Temp) > res.2 <- lm(ozone ~ Solar.R+Wind) > lambda <- -2*(logLik(res.2)-logLik(res.3)) > print(1-pchisq(lambda,1)) log Lik e-09 (df=4) > detach(airquality) lm() 2 glm() nls() 11.8 AIC: AIC L n AIC = 2 ln L + 2n AIC R AIC() extractaic() 2 Akaike s An Information Criterion The (generalized) Akaike *A*n *I*nformation *C*riterion for a fitted parametric model extractaic() MASS S4 step() 40

41 2 AIC AIC(res.3) -2*logLik(res.3)+2*attr(logLik(res.3),"df") extractaic(res.3) res.2 AIC(res.2) extractaic(res.2) AIC AIC() extractaic() step() AIC 44 step(lm(ozone~wind+solar.r+temp)) step() (direction="forward") (direction="backward") (direction="both") direction="backward" step() step() direction ress <- step(lm(ozone~solar.r+wind+temp)) 3 AIC 682 ress AIC step() extractaic() AIC AIC(ress) lm() 111 Ozone Solar.R AIC(ress)-111*(1+log(2*pi)) step() > data(airquality) > attach(airquality) > res <- lm(ozone~solar.r+wind+temp) > ress <- step(res) > summary(ress) > AIC(ress) 5% n n R 8 AIC() 2 ln L + 2θ L θ n p σ n(1 + ln(2πσ 2 )) + 2(p + 1) extractaic() n ln(σ 2 ) + 2p n(1 + ln(2π)) R-help S-plus step.glm() R 41

42 n 1 M.G R step() predict() > data(airquality) > attach(airquality) > res<-lm(ozone~solar.r+wind+temp) > predict(res, list(solar.r=mean(res$model$solar.r), + Wind=mean(res$model$Wind), Temp=mean(res$model$Temp))) [1] > detach(airquality) Solar.R Temp Solar.R mean(solar.r,na.rm=t) Wind= % AIC Wind Ozone Wind Ozone Wind Solar.R 2 > data(airquality) > attach(airquality) > resmr <- nls(ozone ~ a*exp(-b*wind) + c*solar.r, start=list(a=200,b=0.2,c=1)) > summary(resmr) Formula: Ozone ~ a * exp(-b * Wind) + c * Solar.R Parameters: Estimate Std. Error t value Pr(> t ) a e-09 *** b e-11 *** c e-05 *** --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 108 degrees of freedom > AIC(resmr) [1] R 42

43 AIC 2 extractaic() Temp > SRM <- mean(subset(solar.r,!is.na(ozone)&!is.na(solar.r)&!is.na(wind)&!is.na(temp))) > predict(resmr,list(wind=25,solar.r=srm)) ppb Y = β 0 + β 1 X 1 + β 2 X 2 + β 12 X 1 X 2 + ε X 1 Y Y X 2 X 2 X 2 Y (slope) X 1 X 2 Y adjusted mean; X 1 R X 1 C C factor X 2 X Y Y summary(lm(y~c+x)) X C Y C 1 2 C2 (slope) summary(lm(y[c==1]~x[c==1]); summary(lm(y[c==2]~x[c==2]) summary(lm(y~c+x+c:x)) summary(lm(y~c*x)) C X Y Coefficients C2:X R ToothGrowth 10 3 C len supp dose coplot(len~dose supp) 2 43

44 > data(toothgrowth) > attach(toothgrowth) > plot(dose,len,pch=as.integer(supp),ylim=c(0,35)) > legend(max(dose)-0.5,min(len)+1,levels(supp),pch=c(1,2)) > abline(lm1 <- lm(len[supp== VC ]~dose[supp== VC ])) > abline(lm2 <- lm(len[supp== OJ ]~dose[supp== OJ ]),lty=2) > summary(lm1) > summary(lm2) summary(lm1) summary(lm2) > lm3 <- lm(len ~ supp*dose) > summary(lm3) Call: lm(formula = len ~ supp * dose) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) e-09 *** suppvc *** dose e-08 *** suppvc:dose * --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 56 degrees of freedom Multiple R-Squared: , Adjusted R-squared: F-statistic: on 3 and 56 DF, p-value: 6.521e-16 suppvc:dose len % glm() P ln(p/(1 P )) = b 0 + b 1 X b k X k X 1 X 2,...X k X 1 = 0 X 1 = 1 b1 = ln(p 1 /(1 P 1 )) ln(p 0 /(1 P 0 )) = ln(p 1 (1 P 0 )/(P 0 (1 P 1 ))) 44

45 b 1 95% exp(b 1 ± 1.96 SE(b 1 )) library(mass) data(birthwt) Springfield Baystate 189 str(birthwt) low age lwt race smoke ptl ht ui ftv bwt 2.5 kg 1 a (g) a lb. 1 lb kg > require(mass) > data(birthwt) > attach(birthwt) > low <- factor(low) > race <- factor(race, labels=c("white","black","other")) > ptd <- factor(ptl>0) > smoke <- (smoke>0) > ht <- (ht>0) > ui <- (ui>0) > ftv <- factor(ftv) > levels(ftv)[-(1:2)] <- "2+" > bw <- data.frame(low,age,lwt,race,smoke,ptd,ht,ui,ftv) > detach(birthwt) > summary(res <- glm(low ~., family=binomial, data=bw)) > summary(res2 <- step(res)) smoketrue SE % 45

46 > exp( ) [1] > exp( qnorm(0.975)* ) [1] > exp( qnorm(0.975)* ) [1] % [1.08, 5.26] 95% 1 5% 12 survival (1995) Gehan R MASS 2 library() require() 12.1 (1995) p Gehan 6-MP R Surv() summary() summary 46

47 > require(mass) > require(survival) > print(res<-survfit(surv(time,cens)~treat,data=gehan)) Call: survfit(formula = Surv(time, cens) ~ treat, data = gehan) n events median 0.95LCL 0.95UCL treat=6-mp Inf treat=control > summary(res) Call: survfit(formula = Surv(time, cens) ~ treat, data = gehan) treat=6-mp time n.risk n.event survival std.err lower 95% CI upper 95% CI treat=control time n.risk n.event survival std.err lower 95% CI upper 95% CI NA NA NA > plot(res,lty=c(1,2)) > legend(30,0.2,lty=c(1,2),legend=levels(gehan$treat)) LaTeX LaTeX graphicx eps eps OpenOffice.org Draw eps 47

48 difftime() ISOdate() x names dob dod difftime() 4 [x$names=="robert"] Robert alivedays ISOdate(,, ) > x <- data.frame( + names = c("edward","shibasaburo","robert","hideyo"), + dob = c(" "," "," "," "), + dod = c(" "," "," "," ")) > alivedays <- difftime(x$dod,x$dob)[x$names=="robert"] > alivedays/ > difftime(isodate(2005,1,31),x$dob) A 2 B 1 4,6,8,9 2 5,7,12,

49 1 2 i j e ij i d i i j n ij i n i e ij = d i n ij /n i 47 e 11 = 1 n 11 /n 1 = 4/8 = 0.5 i j d ij w i i j u ij u ij = w i (d ij e ij ) 1 u 1 = i (d i1 e i1 ) u 1 = (1 4/8) + (0 3/7) + (1 3/6) + (0 2/5) + (1 2/4) + (1 1/3) + (0 0/2) + (0 0/1) V = V jj = i (n i n ij )n ij d i (n i d i ) n i2 (n i 1) V = (8 4) 4 (7 3) 3 (6 3) 3 (5 2) 2 (4 2) 2 (3 1) *4/64+4*3/49+3*3/36+3*2/25+2*2/16+2*1/ χ 2 = /1.457 = % % R time event group survdiff(surv(time,event)~group) > require(survival) > time2 <- c(4,6,8,9,5,7,12,14) > event <- c(1,1,1,1,1,1,1,1) > group <- c(1,1,1,1,2,2,2,2) > survdiff(surv(time2,event)~group) χ 2 = p = %

50 12.3 z i = (z i1, z i2,..., z ip ) i t h(z i, t) h(z i, t) = h 0 (t) exp(β 1 z i1 + β 2 z i β p z ip ) h 0 (t) t β 1, β 2,..., β p exp(β x z ix ) Cox z i 1 2 t h 0 (t) exp(β 1 z 11 + β 2 z β p z 1p ) exp(β 1 z 21 + β 2 z β p z 2p ) T S(t) T t S(0) = 1 h(t) t Pr(t T < t + t T t) S(t) S(t + t) h(t) = lim = lim t 0 t t 0 ts(t) = ds(t) dt 1 S(t) = d(log(s(t)) dt H(t) = t h(u)du = log S(t) S(t) = exp( H(t)) 0 z S(z, t) H(z, H(z, t) = t 0 h(z, u)du = t 0 h 0 (u) exp(βz)du = exp(βz)h 0 (t) S(z, t) = exp( H(z, t)) = exp{ exp(βz)h 0 (t)} log( log S(z, t)) = βz + log H 0 (t) βz 12.4 β t i t 1 t i 50

51 i L L β L Cox Exact Breslow Efron Exact Breslow R coxph() Efron Breslow Efron Exact time event group coxph(surv(time,event)~group) 2 Exact coxph(surv(time,event)~group, method="exact") Gehan 6-MP > require(mass) > require(survival) > res <- coxph(surv(time,cens)~treat,data=gehan) > summary(res) Call: coxph(formula = Surv(time, cens) ~ treat, data = gehan) n= 42 coef exp(coef) se(coef) z p treatcontrol exp(coef) exp(-coef) lower.95 upper.95 treatcontrol Rsquare= (max possible= ) Likelihood ratio test= 16.4 on 1 df, p=5.26e-05 Wald test = 14.5 on 1 df, p= Score (logrank) test = 17.3 on 1 df, p=3.28e-05 > plot(survfit(res)) 5% 6-MP exp(coef) MP % [2.15, 10.8] 6-MP plot() 2 95% coxph() subset=(x=="6-mp") 51

52 R coxph() strata() time event treat stage coxph(surv(time,event)~treat+strata(stage)) 2 AIC Breslow R Exact Efron 12.6 survreg() R RjpWiki (1995) SAS 2. (2003) R 3. (2004) R 4. (2004) The R Book R 5. (2005) The R Tips R 6. (1995) 7. Armitage P, Berry G, Matthews JNS (2002) Statistical Methods in Medical Research, 4th ed. Blackwell Publishing 8. Maindonald J, Braun J (2003) Data analysis and graphics using R Cambridge Univ. Press 52