# R program for additional exercise 10.3 (VHM 802)

labs <- read.csv("h:/vhm/vhm802/data_csv/hs10_3.csv")
labs$matlab <- labs$lab+labs$mat*0.1

library(nlme)
labs.grp <- groupedData( conc ~ material|lab/matlab, data=labs)
labs.mixed <- lme( conc ~ as.factor(material), random = ~1|lab/matlab, data=labs.grp)
summary(labs.mixed)
anova(labs.mixed)
plot(labs.mixed) # lowest-level residuals
plot(labs.mixed, form=resid(., type="p") ~ fitted(.)|lab, abline=0)
qqnorm(ranef(labs.mixed)$lab[,]) # quantile plot for lab random effects
qqline(ranef(labs.mixed)$lab[,])
qqnorm(ranef(labs.mixed)$matlab[,]) # quantile plot for mat*lab random effects
qqline(ranef(labs.mixed)$matlab[,])

library(lme4)
labs.mixed2 <- lmer( conc ~ as.factor(material) + (1|lab) + (1|matlab), data=labs)
summary(labs.mixed2)