################################################### ### chunk number 1: ################################################### #line 8 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" pp <- read.delim("http://rem.ph.ucla.edu/rob/mld/data/tabdelimiteddata/pain.txt") pp$trt <- as.character(pp$treatment) pp$trt[pp$trial<=3] <- "baseline" pp$trt <- factor(pp$trt, levels=c("baseline","attend","distract","no directions")) pp <- subset(pp, !is.na(pp$l2paintol)) ################################################### ### chunk number 2: ################################################### #line 16 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" source("printnotsoshort.R") library(nlme) library(multcomp) library(lattice) library(latticeExtra) library(reshape) ################################################### ### chunk number 3: ################################################### #line 26 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" p1 <- xyplot(l2paintol~trial | treatment*cs, data=pp, group=id, type="l", col="gray") p1 <- useOuterStrips(p1) plot(p1) ################################################### ### chunk number 4: ################################################### #line 33 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" m <- lme(l2paintol ~ trt*cs, random=~1|id, data=pp) summary(m) ################################################### ### chunk number 5: ################################################### #line 40 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" ids <- unique(pp$id) cook <- numeric(length(ids)) names(cook) <- ids for(i in seq_along(ids)) { mi <- update(m, subset=id!=ids[i]) d <- fixef(mi) - fixef(m) cook[i] <- (d %*% solve(vcov(m)) %*% d)/length(d) } head(sort(cook, decreasing=TRUE)) ################################################### ### chunk number 6: ################################################### #line 52 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" plot(cook) ################################################### ### chunk number 7: ################################################### #line 57 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" red <- c(62, 17, 61, 31, 21, 32) ppX <- pp; ppX$l2paintol[!ppX$id %in% red] <- NA p2 <- xyplot(l2paintol~trial | treatment*cs, data=ppX, group=id, type="l", col="red") plot(p1+p2) ################################################### ### chunk number 8: ################################################### #line 66 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" pp$res <- residuals(m, level=0) ppr <- cast(pp, id ~ trial, value="res")[,-1] ppr <- subset(ppr, complete.cases(ppr)) pc <- princomp(ppr)$scores pc <- as.data.frame(pc) pc$id <- rownames(ppr) head(pc) ################################################### ### chunk number 9: ################################################### #line 76 "/HOME/other/arendahl/courses/EPSY8282-2011/class/class28.Rnw" pcm <- melt(pc, id.vars=5, measure.vars=1:4) head(pcm) p <- xyplot(value~variable, group=id, data=pcm, type="l") plot(p)