# R program for GO Problem 16.4, VHM 802
exams <- read.csv("h:/vhm/vhm802/data_csv/ch16pr4.csv")

# R commands to reproduce Stata listing in Exam 2008, Question 1
library(doBy) # package to facilitate processing by groups
summaryBy( score~ text+exam, data=exams, FUN=mean)
summaryBy( score~ univ+exam, data=exams, FUN=mean)
summaryBy( score~ text+univ, data=exams, FUN=mean)
summaryBy( score~ text, data=exams, FUN=mean)
summaryBy( score~ exam, data=exams, FUN=mean)
summaryBy( score~ univ, data=exams, FUN=mean)
summaryBy( score~ 1, data=exams, FUN=mean)
anova(lm( score~ exam*text+univ, data=exams))
summary (aov( score~ exam*text+univ, data=exams)) # same as above
summary (aov( score~ exam*text+univ+Error(exam:univ), data=exams))
summary (aov( score~ exam*text+univ+Error(text:univ), data=exams))
# note: last two runs produce a warning but give the desired results

# diagnostics for split-plot model analysis by ANOVA-based method
exams.ANOVA <- lm( score~ exam*text+univ+exam:univ, data=exams)
exams.ANOVA.diag <- exams
exams.ANOVA.diag$fit <- fitted(exams.ANOVA)
exams.ANOVA.diag$stdres <- rstandard(exams.ANOVA)
plot(stdres~fit, data=exams.ANOVA.diag)
qqnorm(exams.ANOVA.diag$stdres)
qqline(exams.ANOVA.diag$stdres)
shapiro.test(exams.ANOVA.diag$stdres)

# whole-plot analysis by averaging
exams.means <- summaryBy( score~ univ+exam, data=exams, FUN=mean)
exams.means.ANOVA <- lm( score.mean ~ univ+exam, data=exams.means)
anova(exams.means.ANOVA)
exams.means.diag <- exams.means
exams.means.diag$fit <- fitted(exams.means.ANOVA)
exams.means.diag$stdres <- rstandard(exams.means.ANOVA)
plot(stdres~fit, data=exams.means.diag)
qqnorm(exams.means.diag$stdres)
qqline(exams.means.diag$stdres)
shapiro.test(exams.means.diag$stdres)
# pairwise comparisons between univ based on exams.means.ANOVA omitted here
# see e.g. solution for additional exercise 7.3

# likelihood-based analysis (recommended approach in R)
library(nlme)
exams$semester <- as.numeric(exams$exam)*10+as.numeric(exams$univ)
head(exams) # to check calculation
exams.grp <- groupedData( score ~ exam|semester, data=exams)
exams.mixed <- lme( score ~ exam*text+univ, random = ~1|semester, data=exams)
summary(exams.mixed)
anova(exams.mixed)
plot(exams.mixed) # lowest-level residuals
qqnorm(exams.mixed)
qqnorm(residuals(exams.mixed)) # alternative coding, with best line
qqline(residuals(exams.mixed))
plot(exams.mixed, form=resid(., type="p") ~ fitted(.)|semester, abline=0)
qqnorm(exams.mixed, ~ranef(.)) # quantile plot for semester random effects
qqnorm(ranef(exams.mixed)[,]) # alternative coding, with best line
qqline(ranef(exams.mixed)[,])

# alternative analysis using lme4 library, model specification only
library(lme4)
exams.mixed2 <- lmer( score ~ exam*text+univ + (1|semester), data=exams)
summary(exams.mixed2)