# STAT_516_Lec_09_inclass

ntr <- c(40.89,37.99,37.18,34.98,34.89,42.07,
         41.22,49.42,45.85,50.15,41.99,46.69,
         44.57,52.68,37.61,36.94,46.65,40.23,
         41.90,39.20,43.29,40.45,42.91,39.97)
block <- as.factor(c(1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4))
trt <- as.factor(c(2,5,4,1,6,3,1,3,4,6,5,2,6,3,5,1,2,4,2,4,6,5,3,1))

boxplot(ntr ~ trt)
boxplot(ntr ~ block)

# RCBD : don't include the interaction because
# there is no replication, i.e. n = 1.
lm_fert <- lm(ntr ~ trt + block)
anova(lm_fert) # it's okay, it's a balanced design

# to estimate variance components, use package lmerTest
library(lmerTest)
lmer_fert <- lmer(ntr ~ trt + (1|block))
lmer_fert

ls_means(lmer_fert)
ls_means(lmer_fert,pairwise = TRUE)

qtukey(0.95,6,15)

y4bar <- mean(ntr[trt == 4])
y5bar <- mean(ntr[trt == 5])
me <- qtukey(0.95,6,15) * 2.683 * sqrt(1/4)
lo <- y4bar - y5bar - me
up <- y4bar - y5bar + me
c(lo,up)


# check assumptions of RCBD

plot(lmer_fert)

# Try ignoring the block information.
# then we have just an ordinary one-way ANOVA:

lm_noblock <- lm(ntr ~ trt) # ignore the block
anova(lm_noblock)


skin <- data.frame(resp = c(3,7,9,4,1,5,11,13,8,3,9,12,14,11,5,6,11,12,7,4,5,
                            10,10,6,3,6,12,15,9,5,18,18,15,13,9,7,15,14,9,7),
                   noise = as.factor(c(rep(40,20),rep(80,20))),
                   shock = as.factor(rep(c(rep(.25,5),rep(.5,5),
                                           rep(.75,5),rep(1,5)),2)),
                   subj = as.factor(rep(1:5,8)))
head(skin)

# RCBD with a two-way factorial structure
lm_skin <- lm(resp ~ noise + shock + noise:shock + subj, data = skin)
anova(lm_skin)