# R program for additional exercise 10.3 (VHM 802)

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

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)