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))

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")

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);

reshape

# wide to long
d1<-data.frame(subject=c("an1","an2"),eat1=1:2,eat3=5:6,trt=c("t1","t2"))
d1.1<-reshape(d1, idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",direction="long")
d1.2<-reshape(d1, idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",times=c(1,3),timevar="week",direction="long")

d2<-data.frame(subject=c("an1","an2"),eat1=1:2,eat3=5:6,gain1=14:15,gain3=20:21,trt=c("t1","t2"))
d2.1<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3"),c("gain1","gain3")),v.names=c("eat","gain"),times=c(1,3),timevar="week",direction="long")
d2.2<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",times=c(1,3),timevar="week",direction="long")
d2.3<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3")),drop=c("gain1","gain3"),v.names="eat",times=c(1,3),timevar="week",direction="long")

# long to wide
reshape(d1.2,idvar="subject",timevar="week",v.names="eat",direction="wide")
reshape(d2.1,idvar="subject",timevar="week",v.names=c("eat","gain"),direction="wide")

from Reshaping data with the reshape function in R by Soren Hojsgaard

Paired t-test with Mixed Models

library(MASS)
library(nlme)

# simulated a dataset
paired<-mvrnorm(n=30,mu=c(10,10.2),Sigma=matrix(c(5,4,4,5),2,2))

before<-paired[,1]
after<-paired[,2]


# paired t-test
t.test(before,after,paired=TRUE)



### Mixed Models for a paired t-test 
# reshape data for mixed models
y<-c(before,after)
id<-rep(c(1:length(before)),2)
time<-c(rep(0,length(before)),rep(1,length(after)))


# fit mixed models
summary(lme(y~time, random=~1|id))

Mixed Models - longitudinal data

### Sitka (79 trees with 5 time points)

# random intercept models = compound symmetry
sitka.lme <- lme(size~treat*ordered(Time), random=~1|tree, data=Sitka)

Mixed Models - random coefficient models

library(nlme)

coplot(Y ~ EP|No, data=petrol)

Petrol <- petrol
Petrol[,2:5] <- scale(Petrol[,2:5],scale=F)

# Y~EP for each No; no overal mean
pet1.lm <- lm(Y~No/EP-1, data=Petrol)

# Y~ mu_i + beta*EP; no overal mean
pet2.lm <- lm(Y~No-1+EP, data=Petrol)
pet3.lm <- lm(Y~SG+VP+V10+EP, data=Petrol)

# Random intercept model
pet3.lme <- lme(Y~SG+VP+V10+EP, random=~1|No, data=Petrol)

# To compare likelihood, use "ML"
pet3.lme <- update(pet3.lme, method="ML")
anova(pet3.lme, pet3.lm)
pet4.lme <- update(pet3.lme, fixed=Y~V10+EP)
# pet4.lme <- lme(Y~V10+EP, random=~1|No, method="ML", data=Petrol)
anova(pet3.lme, pet4.lme)

## Random intercept + Random coefficient
pet5.lme <- lme(Y~V10+EP, random=~1+EP|No, method="ML", data=Petrol)
anova(pet4.lme, pet5.lme)