### R code from vignette source 'phia.Rnw' ################################################### ### code chunk number 1: phia.Rnw:13-20 ################################################### options(useFancyQuotes = FALSE, width=85) car.citation <- citation("car") car.citation$key <- "Fox-CAR2011" nlme.citation <- citation("nlme")[1] nlme.citation$key <- "Pinheiro-nlme" lme4.citation <- citation("lme4")[1] lme4.citation$key <- "Bates-lme4" ################################################### ### code chunk number 2: phia.Rnw:24-27 ################################################### toBibtex(car.citation) toBibtex(nlme.citation) toBibtex(lme4.citation) ################################################### ### code chunk number 3: phia.Rnw:238-240 ################################################### library(phia) some(Boik) ################################################### ### code chunk number 4: modboikresid ################################################### mod.boik <- lm(edr ~ therapy*medication, data=Boik) par(mfcol=c(1,2)) plot(mod.boik, 1:2) # Plot diagnostics for the model Anova(mod.boik) ################################################### ### code chunk number 5: phia.Rnw:259-260 ################################################### mod.boik <- lm(edr ~ therapy*medication, data=Boik) par(mfcol=c(1,2)) plot(mod.boik, 1:2) # Plot diagnostics for the model Anova(mod.boik) ################################################### ### code chunk number 6: phia.Rnw:295-296 ################################################### (boik.means <- interactionMeans(mod.boik)) ################################################### ### code chunk number 7: phia.Rnw:309-310 ################################################### interactionMeans(mod.boik, factors="therapy") ################################################### ### code chunk number 8: phia.Rnw:329-330 ################################################### plot(boik.means) ################################################### ### code chunk number 9: phia.Rnw:345-346 ################################################### plot(boik.means, multiple=FALSE) # Not printed in this paper ################################################### ### code chunk number 10: phia.Rnw:416-417 ################################################### testInteractions(mod.boik, fixed="therapy", across="medication") ################################################### ### code chunk number 11: phia.Rnw:473-476 ################################################### boik.mtable <- xtabs(boik.means$"adjusted mean" ~ therapy+medication, boik.means) boik.mtable <- addmargins(boik.mtable, FUN=mean, quiet=TRUE) print(boik.mtable, digits=4) ################################################### ### code chunk number 12: phia.Rnw:484-488 ################################################### boik.resid <- boik.mtable - boik.mtable[4,5] # Subtract the mean boik.resid <- sweep(boik.resid, 1, boik.resid[,5]) # Subtract row means boik.resid <- sweep(boik.resid, 2, boik.resid[4,]) # Subtract column means print(boik.resid, digits=4) ################################################### ### code chunk number 13: phia.Rnw:495-496 ################################################### testInteractions(mod.boik,residual=c("therapy","medication")) ################################################### ### code chunk number 14: residualsplot ################################################### matplot(t(boik.resid[-4,-5]), type="b", xaxt="n", ylab="Interaction residuals") axis(1, at=1:4, labels=levels(Boik$medication)) ################################################### ### code chunk number 15: phia.Rnw:512-513 ################################################### matplot(t(boik.resid[-4,-5]), type="b", xaxt="n", ylab="Interaction residuals") axis(1, at=1:4, labels=levels(Boik$medication)) ################################################### ### code chunk number 16: phia.Rnw:580-581 ################################################### testInteractions(mod.boik, pairwise="therapy", across="medication") ################################################### ### code chunk number 17: phia.Rnw:598-599 ################################################### testInteractions(mod.boik) ################################################### ### code chunk number 18: phia.Rnw:635-642 ################################################### (custom.contr <- contrastCoefficients( therapy ~ control - (T1 + T2)/2, # Control vs. both therapies therapy ~ T1 - T2, # Therapy T1 vs. T2 medication ~ placebo - (D1 + D2 + D3)/3, # Placebo vs. all doses medication ~ D1 - D3, # Min. dose vs. max medication ~ D2 - (D1 + D2 + D3)/3, # Med. dose vs. average data=Boik, normalize=TRUE)) # Normalize to homogeinize the scale ################################################### ### code chunk number 19: phia.Rnw:649-652 ################################################### colnames(custom.contr$therapy) <- c("cntrl.vs.all", "T1.vs.T2") colnames(custom.contr$medication) <- c("plcb.vs.all", "D1.vs.D3", "D2.vs.avg") testInteractions(mod.boik,custom=custom.contr) ################################################### ### code chunk number 20: phia.Rnw:690-691 ################################################### some(OBrienKaiser, 6) ################################################### ### code chunk number 21: phia.Rnw:703-704 ################################################### (idata <- expand.grid(hour=ordered(1:5), phase=c("pre", "post","fup"))) ################################################### ### code chunk number 22: phia.Rnw:710-711 ################################################### addmargins(table(OBrienKaiser[c("gender","treatment")])) ################################################### ### code chunk number 23: phia.Rnw:717-721 ################################################### mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment*gender, data=OBrienKaiser) ################################################### ### code chunk number 24: phia.Rnw:730-731 ################################################### Anova(mod.ok, idata=idata, idesign=~phase*hour, type=3) ################################################### ### code chunk number 25: okplot ################################################### ok.means <- interactionMeans(mod.ok, c("treatment","phase"), idata=idata) plot(ok.means, errorbar="ci95") ################################################### ### code chunk number 26: phia.Rnw:766-767 ################################################### ok.means <- interactionMeans(mod.ok, c("treatment","phase"), idata=idata) plot(ok.means, errorbar="ci95") ################################################### ### code chunk number 27: phia.Rnw:785-786 ################################################### testInteractions(mod.ok, pairwise=c("treatment", "phase"), idata=idata) ################################################### ### code chunk number 28: phia.Rnw:861-862 ################################################### str(Prestige) ################################################### ### code chunk number 29: phia.Rnw:875-877 ################################################### mod.prestige <- lm(prestige ~ (log2(income)+education)*type, data=Prestige) Anova(mod.prestige) ################################################### ### code chunk number 30: prestigeplot ################################################### plot(interactionMeans(mod.prestige, atx="type", slope="log2(income)")) testInteractions(mod.prestige, pairwise="type", slope="log2(income)") ################################################### ### code chunk number 31: phia.Rnw:896-897 ################################################### plot(interactionMeans(mod.prestige, atx="type", slope="log2(income)")) testInteractions(mod.prestige, pairwise="type", slope="log2(income)") ################################################### ### code chunk number 32: phia.Rnw:915-916 ################################################### table(Prestige$type) # Frequencies of occupation types ################################################### ### code chunk number 33: phia.Rnw:940-943 ################################################### mod.prestige2 <- update(mod.prestige, formula=.~.+(log2(income):education)*type, subset=(Prestige$type!="wc")) Anova(mod.prestige2) ################################################### ### code chunk number 34: phia.Rnw:951-952 ################################################### testInteractions(mod.prestige2, pairwise="type", slope="log2(income)") ################################################### ### code chunk number 35: phia.Rnw:966-970 ################################################### # Look quantiles of the model frame (a subset of the original data) quantile(model.frame(mod.prestige2)$education) testInteractions(mod.prestige2, pairwise="type", slope="log2(income)", covariates=c(education=14)) ################################################### ### code chunk number 36: phia.Rnw:1031-1032 ################################################### str(AMSsurvey) ################################################### ### code chunk number 37: phia.Rnw:1055-1056 ################################################### ftable(xtabs(count ~ sex + citizen + type, AMSsurvey)) ################################################### ### code chunk number 38: phia.Rnw:1066-1067 ################################################### mod.ams <- glm(count ~ type*(sex+citizen), family=poisson, data=AMSsurvey) ################################################### ### code chunk number 39: phia.Rnw:1073-1074 ################################################### Anova(mod.ams) ################################################### ### code chunk number 40: amsplot ################################################### ams.means <- interactionMeans(mod.ams) plot(ams.means, atx="type", traces=c("sex","citizen")) ################################################### ### code chunk number 41: phia.Rnw:1092-1093 ################################################### ams.means <- interactionMeans(mod.ams) plot(ams.means, atx="type", traces=c("sex","citizen")) ################################################### ### code chunk number 42: phia.Rnw:1125-1127 ################################################### testInteractions(mod.ams, pairwise=c("type","sex")) # test type:sex testInteractions(mod.ams, pairwise=c("type","citizen")) #test type:citizen ################################################### ### code chunk number 43: phia.Rnw:1176-1177 (eval = FALSE) ################################################### ## dm.de <- family(model)$mu.eta ################################################### ### code chunk number 44: phia.Rnw:1302-1308 ################################################### Snijders <- nlme::bdf[c("langPRET","langPOST", # Outcomes "pupilNR", "IQ.ver.cen", "ses", "sex", # Student-related variables "schoolNR", "schoolSES", "avg.IQ.ver.cen")] # School-related variables Snijders$sex <- factor(Snijders$sex, labels=c("F","M")) names(Snijders) <- c("score.1","score.2","student","IQ","SES","sex","school","avgSES","avgIQ") ################################################### ### code chunk number 45: phia.Rnw:1315-1317 ################################################### Snijders$IQ2.pos <- with(Snijders, (IQ > 0)*IQ^2) Snijders$IQ2.neg <- with(Snijders, (IQ < 0)*IQ^2) ################################################### ### code chunk number 46: phia.Rnw:1328-1329 ################################################### library(lme4) ################################################### ### code chunk number 47: phia.Rnw:1331-1335 (eval = FALSE) ################################################### ## form1 <- score.2 ~ ## IQ * SES + IQ2.pos + IQ2.neg + sex + avgIQ * avgSES + # Fixed part ## (IQ | school) # Random part ## mod.snijders.1 <- lmer(form1, data=Snijders) ################################################### ### code chunk number 48: phia.Rnw:1347-1349 (eval = FALSE) ################################################### ## form2.1 <- update(form1, (score.1+score.2)/2~.) ## form2.2 <- update(form1, (score.2-score.1)~.) ################################################### ### code chunk number 49: phia.Rnw:1368-1374 ################################################### Snijders.long <- reshape(Snijders, direction="long", idvar="student", varying=list(c("score.1","score.2")), v.names="score", timevar="repetition") # The within-subjects factor must be coded as a factor Snijders.long$repetition <- as.factor(Snijders.long$repetition) # See the variables of the long data frame str(Snijders.long) ################################################### ### code chunk number 50: phia.Rnw:1390-1397 ################################################### form3 <- score ~ repetition * (IQ * SES + IQ2.pos + IQ2.neg + sex) + # Student-related avgIQ * avgSES + # School-related (IQ | school) + (1 | student) # Random part mod.snijders.3 <- lmer(form3, data=Snijders.long) # See the main parameters of the model (ommit correlations table) print(mod.snijders.3, correlation=FALSE) ################################################### ### code chunk number 51: phia.Rnw:1428-1429 ################################################### Anova(mod.snijders.3) ################################################### ### code chunk number 52: phia.Rnw:1442-1448 ################################################### # Cell means interactionMeans(mod.snijders.3) # Simple effect of sex at each repetition testInteractions(mod.snijders.3, fixed="repetition", across="sex") # Pairwise interactions (default) testInteractions(mod.snijders.3) ################################################### ### code chunk number 53: phia.Rnw:1481-1482 ################################################### graphics.off()