############################Split-plot experiment########################## library(lmerTest) #this is the package suggested library(nlme) #for random effects, degrees of freedom has some problem require(ggplot2) require(reshape2) #to change long format to wide format or change from wide format to long format library(heplots) #test equivalence of covariance matrix library(multcomp) #for multiple comparison library(xtable) library(lsmeans) library(multcompView) #hrate is in long format hrate<-read.table("C:/jenn/teaching/stat579/data/heartrate.txt",header=TRUE) nrow(hrate) #120 hrate[1:20,] hrate$time<-factor(hrate$time) hrate$drug<-factor(hrate$drug) hrate$subject<-factor(hrate$subject) boxplot(rate~drug, data=hrate, main="boxplot of rate vs drug") #use ggplot for rate vs drug for each time period p <- ggplot(hrate, aes(factor(time), rate)) + geom_boxplot(aes(fill = drug)) p ## d2 is a data set in wide format, with 30 variables and 6 columns, drug, #subject, time1-time4, you can change d2 to the long format, such as the data #hrate with 120 obsn and 4 columns drug, subject, time, and rate d2<-read.table("C:/jenn/teaching/stat579/data/splitd2.txt",header=TRUE) d2.long <- melt(d2, # id.vars: ID variables # all variables to keep but not split apart on id.vars=c("subject","drug"), # measure.vars: The source columns # (if unspecified then all other variables are measure.vars) measure.vars = c("time1", "time2", "time3", "time4"), # variable.name: Name of the destination column identifying each # original column that the measurement came from variable.name = "time", # value.name: column name for values in table value.name = "rate" ) nrow(d2) d2[1:10,] nrow(d2.long) d2.long[1:20,] d2p<-d2[1:10,3:6] d2a<-d2[11:20,3:6] d2b<-d2[21:30,3:6] cov(d2p) cov(d2a) cov(d2b) ##test if the three covariance matrices are equivalent across the groups tres <- boxM(d2[, 3:6], d2[, "drug"]) tres summary(tres) ##checking interactions #seems interaction is significant interaction.plot(hrate$drug,hrate$time,hrate$rate) #ANOVA table lmod <- aov(rate ~ drug*time + Error(factor(subject):factor(drug)), data=hrate) summary(lmod) ##fit full model with two treatments ##note that there are 30 subjects with 1-30, identifiable betweent the groups ##adding nested structure or not doesn't make difference myfit<-lmer(rate~drug*time+(1|subject),data=hrate) summary(myfit) print(xtable(anova(myfit))) #this is for the latex version to enter the table anova(myfit) myfit2<-lmer(rate~drug*time+(1|subject/drug),data=hrate) summary(myfit2) anova(myfit2) #same result if using d2.long data myfit3<-lmer(rate~drug*time+(1|subject),data=d2.long) summary(myfit3) anova(myfit3) #however, if within each group p, a and b, subjects are listed as 1-10, 1-10, 1-10 respectively #you will need to add the nested structure, otherwise, R will consider subject 1 from group p, a and b as the #same subject #d3 is constructed as "subjects are listed as 1-10, 1-10, 1-10 respectively within the three groups" d3<-read.table("C:/jenn/teaching/stat579/data/splitd3.txt",header=TRUE) d3.long <- melt(d3, # id.vars: ID variables # all variables to keep but not split apart on id.vars=c("subject","drug"), # measure.vars: The source columns # (if unspecified then all other variables are measure.vars) measure.vars = c("time1", "time2", "time3", "time4"), # variable.name: Name of the destination column identifying each # original column that the measurement came from variable.name = "time", # value.name: column name for values in table value.name = "rate" ) nrow(d3) d3 #in this case, without nested structure, R analyze as a completely randomized block design myfit4<-lmer(rate~drug*time+(1|subject),data=d3.long) summary(myfit4) anova(myfit4) #adding nested structure myfit5<-lmer(rate~drug*time+(1|subject/drug),data=d3.long) #subject nested in drug summary(myfit5) anova(myfit5) ##wrong degrees of freedom from nlme package again, note that the degrees of #freedom for testing drug effect is 18 instead of 27 myfit6<-lme(rate~drug*time,random=~1|subject/drug,data=d3.long) summary(myfit6) anova(myfit6) #We could check of if the random effect is important or not, by the step #function of lmerTest: (usually, people would prefer to include the random effect #in the model whether they are signnificant or not) mystep <- step(myfit) #subject is significant, need to keep mystep$rand.table mystep$anova.table mystep$diffs.lsmeans.table #mutlitple comparisons ##since interaction is significant, let's do multiple comparisons marginal<-lsmeans(myfit,pairwise~drug:time) cld(marginal, alpha=0.05, Letters=letters, adjust="tukey") ##diagnostics fit.lmer<-predict(myfit) resid.lmer<-resid(myfit) par(mfrow=c(2,2)) plot(resid.lmer~fit.lmer,main="fitted values vs. residuals",xlab="fitted",ylab="residuals") abline(0,0,lty=2) qqnorm(resid.lmer) qqline(resid.lmer) hist(resid.lmer,main="histogram of residuals",xlab="residuals") #random effects normality par(mfrow=c(1,2)) alphahat<-ranef(myfit) hist(alphahat[,1]) #you may need to use #hist(alphahat$subject[,1]) abline(0,0,lty=2) qqnorm(alphahat[,1]) qqline(alphahat[,1])