library(lmerTest)

?Pastes ## info about data set
str(Pastes)
head(Pastes, 10)

## Visualize data (code taken from help file) ####
require(lattice)
dotplot(cask ~ strength | reorder(batch, strength), Pastes,
        strip = FALSE, strip.left = TRUE, layout = c(1, 10),
        ylab = "Cask within batch",
        xlab = "Paste strength", jitter.y = TRUE)
## Modifying the factors to enhance the plot
Pastes <- within(Pastes, batch <- reorder(batch, strength))
Pastes <- within(Pastes, sample <- reorder(reorder(sample, strength),
                                           as.numeric(batch)))
dotplot(sample ~ strength | batch, Pastes,
        strip = FALSE, strip.left = TRUE, layout = c(1, 10),
        scales = list(y = list(relation = "free")),
        ylab = "Sample within batch",
        xlab = "Paste strength", jitter.y = TRUE)

## Get overview of number of observations ####
table(Pastes$batch, Pastes$cask)

## classical appraoch with aov####

## model everything as fixed
fit <- aov(strength ~ batch + cask %in% batch, data = Pastes)
#fit <- aov(strength ~ batch/cask, data = Pastes) ## equivalent formulation
summary(fit)

pf(27.489 / 17.545, 9, 20, lower.tail = FALSE)

## lmer approaches (all equivalent) ####
fm1 <- lmer(strength ~ (1 | batch/cask), data = Pastes)
summary(fm1)
confint(fm1, oldNames = FALSE)

fm2 <- lmer(strength ~ (1 | batch) + (1 | cask:batch), data = Pastes)
rand(fm2)

## sample is already encoded in "nested" form
fm3 <- lmer(strength ~ (1 | batch) + (1 | sample), Pastes)

## the following is *not* appropriate here, why?
(fm4 <- lmer(strength ~ (1 | batch) + (1 | cask), Pastes)) 

## residual analysis ####
plot(fm1) ## Tukey-Anscombe plot
par(mfrow = c(2,2))
qqnorm(ranef(fm1)$batch[,1], main = "batch")
qqnorm(ranef(fm1)$'cask:batch'[,1], main = "cask")
qqnorm(resid(fm1), main = "residuals")