# STAT_516_Lec_08_inclass

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)

lm_machines <- lm(y ~ machine)
anova_machines <- anova(lm_machines)

MStrt <- anova_machines$`Mean Sq`[1]
MSerr <- anova_machines$`Mean Sq`[2]

n <- 5
sge2_hat <- MSerr # this is my estimate of sigme_epsilon^2
sgA2_hat <- (MStrt - MSerr)/n

sgA2_hat
sge2_hat

# install.packages("lmerTest")
library(lmerTest)

# lmer function can handle random effects models
lmer_machines <- lmer(y ~ 1 + (1|machine))
summary(lmer_machines)

predict(lmer_machines, random.only=TRUE)

# qq plot of residuals
yhat <- predict(lmer_machines)
ehat <- y - yhat
qqnorm(ehat)

# residuals versus fitted values plot
plot(lmer_machines)


# get the predicted values of the Ai

y.. <- mean(y)
yi. <- aggregate(y ~ machine, FUN=mean)[,2]
Ai_hat <- n*sgA2_hat / (n*sgA2_hat + sge2_hat)* (yi. - y..)
Ai_hat
abline(h = y..)