################################################### ### chunk number 1: ################################################### #line 12 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" library(nlme) library(lattice) library(latticeExtra) library(reshape) library(multcomp) black <- standard.theme("pdf") black$superpose.line$col <- "black" black$superpose.symbol$col <- "black" black$plot.symbol$col <- "black" black$plot.line$col <- "black" options(show.signif.stars=FALSE) source("printshorter2.R") ################################################### ### chunk number 2: ################################################### #line 27 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" ravens <- read.delim("http://rem.ph.ucla.edu/rob/mld/data/tabdelimiteddata/cognitive.txt") r1 <- subset(ravens, select=c("id","rn","relmonth", "sex","ses1","ravens","treatment")) r1 <- r1[complete.cases(r1),] r1$id <- factor(r1$id) r1$sex <- factor(r1$sex, levels=c("girl","boy")) r1$treatment <- factor(r1$treatment, levels=c("control","milk","calorie","meat")) r1$trt <- r1$treatment r1$trt[r1$rn==1] <- "control" ################################################### ### chunk number 3: ################################################### #line 47 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" rmeans <- melt(r1, id.vars=c("id","rn","sex","treatment"), measure.vars=c("relmonth","ravens")) rmeans <- cast(rmeans, rn+treatment~variable, fun.aggregate=mean) black2 <- black black2$superpose.line$lty <- c(1:4) p <- xyplot(ravens~relmonth, group=treatment, type=c("p","l"), data=rmeans, auto.key=list(space="right", points=FALSE, lines=TRUE), par.settings=black2) ################################################### ### chunk number 4: ################################################### #line 55 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" plot(p) ################################################### ### chunk number 5: ################################################### #line 59 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" m1 <- lme(ravens~relmonth+relmonth:trt, random=~1|id, data=r1) summary(m1) ################################################### ### chunk number 6: ################################################### #line 64 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" m0 <- lme(ravens~relmonth, random=~1|id, data=r1) m0b <- update(m0, method="ML") m1b <- update(m1, method="ML") anova(m0b, m1b) ################################################### ### chunk number 7: ################################################### #line 71 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" K1 <- rbind(`control slope`=c(0,1,0,0,0), `milk slope` =c(0,1,1,0,0), `calorie slope`=c(0,1,0,1,0), `meat slope` =c(0,1,0,0,1)) K2 <- rbind(`meat-control`=c(0,0,0,0,1), `meat-calorie`=c(0,0,0,-1,1), `meat-milk`=c(0,0,-1,0,1), `control-calorie`=c(0,0,0,-1,0), `control-milk`=c(0,0,-1,0,0), `calorie-milk`=c(0,0,-1,1,0)) ################################################### ### chunk number 8: ################################################### #line 85 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" K1 t1 <- glht(m1, linfct=K1) summary(t1, test=adjusted(type="none")) ################################################### ### chunk number 9: ################################################### #line 91 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" foo <- expand.grid(relmonth=c(-2,-0.0001,0.0001,25),treatment=levels(r1$trt)) foo$trt <- foo$treatment foo$trt[foo$relmonth<=0] <- "control" foo$ravens <- 1 foo$y0 <- model.matrix(formula(m0), data=foo) %*% fixef(m0) foo$y1 <- model.matrix(formula(m1), data=foo) %*% fixef(m1) p <- xyplot(ravens~relmonth|treatment, type="p", data=rmeans, par.settings=black) p0 <- xyplot(y0~relmonth|treatment, foo, type="l", col="black", lty=2) p1 <- xyplot(y1~relmonth|treatment, foo, type="l", col="black") ################################################### ### chunk number 10: ################################################### #line 105 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" plot(p+p0+p1) ################################################### ### chunk number 11: ################################################### #line 116 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" K2 t2 <- glht(m1, linfct=K2) summary(t2) ################################################### ### chunk number 12: ################################################### #line 126 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" m2 <- lme(ravens~relmonth+trt, random=~1|id, data=r1) summary(m2) ################################################### ### chunk number 13: ################################################### #line 132 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" m3 <- lme(ravens~relmonth*trt, random=~1|id, data=r1) summary(m3) ################################################### ### chunk number 14: ################################################### #line 138 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" #m4 <- lme(ravens~relmonth*treatment, random=~1|id, data=r1) #summary(m4) #m5 <- lme(ravens~relmonth*treatment + relmonth:trt, random=~1|id, data=r1) #summary(m4) ################################################### ### chunk number 15: ################################################### #line 145 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" m2b <- update(m2, method="ML") m3b <- update(m3, method="ML") #m4b <- update(m4, method="ML") #m5b <- update(m5, method="ML") anova(m0b,m1b, m2b,m3b) anova(m1b,m3b) ################################################### ### chunk number 16: ################################################### #line 154 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" foo$y2 <- model.matrix(formula(m2), data=foo) %*% fixef(m2) foo$y3 <- model.matrix(formula(m3), data=foo) %*% fixef(m3) #foo$y4 <- model.matrix(formula(m4), data=foo) %*% fixef(m4) #foo$y5 <- model.matrix(formula(m5), data=foo) %*% fixef(m5) p1 <- xyplot(y1~relmonth, group=treatment, foo, type="l", par.settings=black2) p2 <- xyplot(y2~relmonth, group=treatment, foo, type="l", par.settings=black2) p3 <- xyplot(y3~relmonth, group=treatment, foo, type="l", par.settings=black2) #p4 <- xyplot(y4~relmonth, group=treatment, foo, type="l", par.settings=black2) #p5 <- xyplot(y5~relmonth, group=treatment, foo, type="l", par.settings=black2) ps <- c(m1=p1,m2=p2,m3=p3)#,m4=p4,m5=p5) ps <- update(ps, as.table=TRUE, scale=list(y="same")) ################################################### ### chunk number 17: ################################################### #line 169 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" plot(ps) ################################################### ### chunk number 18: ################################################### #line 177 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" s1 <- lme(ravens~relmonth+relmonth:trt+sex*relmonth, random=~1|id, data=r1) summary(s1) ################################################### ### chunk number 19: ################################################### #line 183 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" s2 <- lme(ravens~relmonth+relmonth:trt+sex*relmonth + relmonth:trt:sex, random=~1|id, data=r1) summary(s2) ################################################### ### chunk number 20: ################################################### #line 190 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" s1b <- update(s1, method="ML") s2b <- update(s2, method="ML") anova(m1b,s1b,s2b) ################################################### ### chunk number 21: ################################################### #line 196 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" foo <- expand.grid(relmonth=c(-2,-0.0001,0.0001,25),treatment=levels(r1$trt), sex=levels(r1$sex)) foo$trt <- foo$treatment foo$trt[foo$relmonth<=0] <- "control" foo$ravens <- 1 rmeans <- melt(r1, id.vars=c("id","rn","sex","treatment"), measure.vars=c("relmonth","ravens")) rmeans <- cast(rmeans, rn+sex+treatment~variable, fun.aggregate=mean) foo$y1 <- model.matrix(formula(m1), data=foo) %*% fixef(m1) foo$yy1 <- model.matrix(formula(s1), data=foo) %*% fixef(s1) foo$yy2 <- model.matrix(formula(s2), data=foo) %*% fixef(s2) p <- useOuterStrips(xyplot(ravens~relmonth|treatment*sex, type=c("p"), data=rmeans, par.settings=black2, lty=1)) p0 <- useOuterStrips(xyplot(y1~relmonth|treatment*sex, foo, type="b", col="gray")) p1 <- useOuterStrips(xyplot(yy1~relmonth|treatment*sex, foo, type="b", col="black")) p2 <- useOuterStrips(xyplot(yy2~relmonth|treatment*sex, foo, type="b", col="black", lty=2)) ################################################### ### chunk number 22: ################################################### #line 216 "/HOME/other/arendahl/courses/EPSY8282-2011/class/case-cognitive2.Rnw" plot(p+p0+p1+p2)