# R program for additional exercise 10.4 (VHM 802)

mice <- read.csv("h:/vhm/vhm802/data_csv/hs10_4.csv")

# anova-based analysis (anova table only, no variance components)
mice.nested <- lm( pH ~ as.factor(strain)/as.factor(litter), data=mice)
summary(mice.nested)
anova(mice.nested)

# likelihood-based analysis
library(nlme)
mice.grp <- groupedData( pH ~ strain|litterid, data=mice)
mice.mixed <- lme( pH ~ as.factor(strain), random = ~1|litterid, data=mice.grp)
summary(mice.mixed)
anova(mice.mixed)
plot(mice.mixed) # lowest-level residuals
plot(mice.mixed, form=resid(., type="p") ~ fitted(.)|litterid, abline=0)
qqnorm(ranef(mice.mixed)) # quantile plot for litterid random effects
qqline(ranef(mice.mixed)[,])

library(lme4)
mice.mixed2 <- lmer( pH ~ as.factor(strain) + (1|litterid), data=mice)
summary(mice.mixed2)