# STAT_516_Lec_09

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))


stripchart(ntr ~ trt, vertical=T, pch = 1)
stripchart(ntr ~ block, vertical=T, pch = 1)


# compute the SS for these data:

a <- 6
b <- 4

y <- ntr
y.. <- predict(lm(y ~ 1))
yi. <- predict(lm(y ~ trt))
y.j <- predict(lm(y ~ block))

SST <- sum((y - y..)**2)
SSA <- sum((yi. - y..)**2)
SSB <- sum((y.j - y..)**2)
SSE <- sum((y - (yi. + y.j - y..))**2)

# these are equal
SSA + SSB + SSE
SST

MSA <- SSA / (a - 1)
MSB <- SSB / (b - 1)
MSE <- SSE / ((a-1)*(b-1))

FA <- MSA / MSE
FB <- MSB / MSE


pA <- 1 - pf(FA,a-1,(a-1)*(b-1))
pB <- 1 - pf(FB,b-1,(a-1)*(b-1))

pA
pB

# can get anova table (correct!)
lm_out <- lm(ntr ~ trt + block)
anova(lm_out)

# this won't work because we have no replications
lm_out_bad <- lm(ntr ~ trt + block + trt:block)
anova(lm_out_bad)


# estimate the variance components
sigma_B_sq <- (MSB - MSE) / a
sigma_e_sq <- MSE

sigma_B <- sqrt(sigma_B_sq)
sigma_e <- sqrt(sigma_e_sq)

sigma_B
sigma_e

#  use lmer() function

library(lmerTest)
lmer_out <- lmer(ntr ~ trt + (1|block))
lmer_out


ls_means(lmer_out)
ls_means(lmer_out,pairwise = T)

aggregate(ntr ~ trt, FUN=mean)

# what if we ignore the blocks??

anova(lm(ntr ~ trt))


# RCBD with replication:

y <- c(104,114,109,124,
       134,130,154,164,
       146,142,152,156,
       147,160,160,163,
       133,146,156,161)

nitr <- as.factor(c(0,0,0,0,
                    50,50,50,50,
                    100,100,100,100,
                    150,150,150,150,
                    200,200,200,200))

blk <- as.factor(c(1,1,2,2,
                   1,1,2,2,
                   1,1,2,2,
                   1,1,2,2,
                   1,1,2,2))

boxplot(y~ nitr)

a <- 5
b <- 2
n <- 2

yi.. <- predict(lm(y ~ nitr))
y.j. <- predict(lm(y ~ blk))
# yij. <- predict(lm(y ~ nitr + blk + nitr:blk))

y... <- mean(y)
SStot <- sum( (y - y...)**2)
SStrt <- sum((yi.. - y...)**2)
SSblk <- sum((y.j. - y...)**2)
SSerror <- sum((y - (yi.. + y.j. - y...))**2)

MStrt <- SStrt / (a-1)
MSblk <- SSblk / (b-1)
MSerror <- SSerror / (a*b*n - a - b + 1)


Ftrt <- MStrt / MSerror
Fblk <- MSblk / MSerror

ptrt <- 1 - pf(Ftrt,a - 1, a*b*n - a - b + 1)
pblk <- 1 - pf(Fblk,b - 1, a*b*n - a - b + 1)


# this table is correct when we do not have the interaction term
anova(lm(y ~ nitr + blk))


# estimate variance components:

sigma_B_sq <- (MSblk - MSerror) / (n*a)
sigma_e_sq <- MSerror

sigma_B <- sqrt(sigma_B_sq)
sigma_e <- sqrt(sigma_e_sq)

# lmer()


lmer_out <- lmer(y ~ nitr + (1|blk))
lmer_out

# confidence interval for difference between two groups:


y50.. <- mean(y[nitr == "50"])
y0.. <- mean(y[nitr == "0"])
alpha <- 0.05
tval <- qt(1 - alpha/2,a*b*n - a -b + 1)
me <- tval * sqrt(MSerror) * sqrt(2/(n*b))
lo <- y50.. - y0.. - me
up <- y50.. - y0.. + me
c(lo,up)

confint(lmer_out)
ls_means(lmer_out,pairwise=T)