# brief R program for additional exercise 10.5 (VHM 802)
# focus is on modelling and presentation

hib <- read.csv("h:/vhm/vhm802/data_csv/hs10_5.csv")
hib$troughid <- 10*as.numeric(hib$water) + hib$trough

# hierarchical model analysis using nlme library
library(nlme)
hib.grp <- groupedData( weight ~ var|trough, outer= ~water, data=hib)
plot(hib.grp, outer = ~ water)
hib.hier <- lme( weight ~ var*water, random=~1|troughid, data=hib.grp)
summary(hib.hier)
anova(hib.hier)
plot(hib.hier) # lowest-level residuals
plot(hib.hier, form=resid(., type="p") ~ fitted(.)|troughid, abline=0)
qqnorm(hib.hier)
qqnorm(residuals(hib.hier))
qqline(residuals(hib.hier))
qqnorm(hib.hier, ~ranef(.)) # quantile plot for trough random effects
qqnorm(ranef(hib.hier)[,]) 
qqline(ranef(hib.hier)[,])

# least squares means with (correct) SE, using add-on libraries
library(estimability)
library(lsmeans)
lsmeans(hib.hier, "water")
lsmeans(hib.hier, "var")
lsmeans(hib.hier, ~ var|water)
# pairwise comparisons with correct SE and z-based P-values, with Tukey adjustments
lsmeans(hib.hier, pairwise ~ water)
lsmeans(hib.hier, pairwise ~ var|water)
lsmeans(hib.hier, pairwise ~ water|var)

# hierarchical model analysis using lme4 library
library(lme4)
hib.hier2 <- lmer( weight ~ var*water + (1|troughid), data=hib.grp)
summary(hib.hier2)