# STAT_516_Lec_08

y <- c(128,127,129,126,128,121,120,123,122,125,
       126,125,127,125,124,125,126,129,128,127)
machine <- as.factor(c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4))

boxplot(y ~ machine)

# get SS as in fixed effects models

lm_out <- lm(y ~ machine)
y.. <- mean(y)
SStot <- sum((y - y..)**2)
SSerror <- sum(lm_out$residuals**2)
SStrt <- SStot - SSerror

a <- 4
n <- 5
MStrt <- SStrt / (a-1)
MSerror <- SSerror / (a*(n-1))

Ftest <- MStrt / MSerror
pval <- 1 - pf(Ftest,a-1, a*(n-1))

anova(lm_out)

# Get method of moments estimators for the two variance components

sigma_A_sq <- (MStrt - MSerror)/n
sigma_e_sq <- MSerror

sigma_A <- sqrt(sigma_A_sq)
sigma_e <- sqrt(sigma_e_sq)

sigma_A
sigma_e

# lm() function is not set up for random effects

library(lmerTest)
# install.packages("lmerTest") # the first time


lmer_out <- lmer(y ~ 1 + (1|machine))
summary(lmer_out)

predict(lmer_out, random.only = T)

# confidence interval for the overall mean

alpha <- 0.05
tval <- qt(1-alpha/2,a - 1)
lo <- y.. - tval * sqrt(MStrt/(a*n))
up <- y.. + tval * sqrt(MStrt/(a*n))
c(lo,up)

# CI for overall mean (slightly different one, but fine)


confint(lmer_out,level = 0.95)


# checking assumptions

# fitted values
yhat <- predict(lmer_out)
ehat <- y - yhat

qqnorm(ehat)

# residuals versus fitted values
plot(ehat ~ yhat)

# automatic residuals versus fitted values plot
plot(lmer_out)