Monday, March 30, 2009

McNemar's test

# simulate a dataset
group<-rep(c('case','cont'),rep(100,2))
pre<-rbinom(200,1,0.7)
post<-rbinom(200,1,0.8)
mcnemar.test(pre,post)

# create a variable post - pre
matched<-data.frame(pre,post,group)
matched$change<-post-pre

# McNemar's test from logistic regression (0 has to be deleted)
summary(glm(factor(change)~1, data=matched[matched$change!=0,],family=binomial))

summary(glm(factor(change)~group, data=matched[matched$change!=0,],family=binomial))

Friday, March 27, 2009

Multiple Comparisons

### Use warpbreaks dataset
amod<-aov(breaks~tension, data=warpbreaks)

### 1-way ANOVA (exact)
library(multcomp)
summary(glht(amod, linfct=mcp(tension="Dunnett")))
summary(glht(amod, linfct=mcp(tension="Tukey")))

# CI
confint(glht(amod, linfct=mcp(tension="Dunnett")), level=0.9)

# create a contrast
contr<-rbind("M-L"=c(-1,1,0), "H-L"=c(-1,0,1), "H-M"=c(0,-1,1))
glht(amod, linfct=mcp(tension=contr))

# Tukey test
TukeyHSD(chickwts.out,"feed")

# Closed test
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="none")
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="bonf")
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="hommel")

Wednesday, March 25, 2009

ROC curve

library(ROC)

R1 <- rocdemo.sca( rbinom(40,1,.3), rnorm(40), caseLabel="new case", markerLabel="demo Marker" );

plot(R1, line=TRUE, show.thresh=TRUE);

plot(R1, line=TRUE);