#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#TWO-WAY ANALYSIS OF VARIANCE
#ONE OBSERVATION PER CELL
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#CONSTRUCTION OF DATA
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 blend_rep(paste("Blend",1:5),4)
 treatment_rep(LETTERS[1:4],rep(5,4))
 yield_scan()
 89 84 81 87 79
 88 77 87 92 81
 97 92 87 89 80
 94 79 85 84 88
 
 pen.df_data.frame(blend,treatment,yield)
 pen.df

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#DIAGNOSTIC PLOTS FOR THE MODEL
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 plot.design(pen.df)
 plot.design(pen.df,fun=median)

 par(mfrow=c(1,2))
 plot.factor(pen.df)

 attach(pen.df)
 interaction.plot(treatment, blend, yield)
 interaction.plot(blend, treatment, yield)
 detach()

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#TABLE OF ANALYSIS OF VARIANCE
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 pen.aov_aov(yield~blend+treatment,pen.df)
 summary(pen.aov)

 model.tables(pen.aov)
 model.tables(pen.aov,type="means")

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#DIAGNOSTIC PLOTS FOR RESIDUALS
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 par(mfrow=c(1,1))
 hist(resid(pen.aov))
 par(mfrow=c(1,2))
 qqnorm(resid(pen.aov))
 plot(fitted(pen.aov),resid(pen.aov))


 tukey.1_function(aov.obj, data)
   {
            vnames <- names(aov.obj$contrasts)
            if(length(vnames) != 2)
                stop("the model must be two-way")
            vara <- data[, vnames[1]]
            varb <- data[, vnames[2]]
            na <- length(levels(vara))
            nb <- length(levels(varb))
            resp <- data[, as.character(attr(aov.obj$terms,
                   "variables")[attr(aov.obj$terms, "response" )])]
            cfs <- coef(aov.obj)
            alpha.A <- aov.obj$contrasts[[vnames[1]]] %*% cfs[
                       aov.obj$assign[[vnames[1]]]]
            alpha.B <- aov.obj$contrasts[[vnames[2]]] %*% cfs[
                       aov.obj$assign[[vnames[2]]]]
            r.mat <- matrix(0, nb, na)
            r.mat[cbind(as.vector(unclass(varb)), as.vector(
                 unclass(vara)))] <- resp
            SS.theta.num <- sum((alpha.B %*% t(alpha.A)) *
                            r.mat)^2
            SS.theta.den <- sum(alpha.A^2) * sum(alpha.B^2)
            SS.theta <- SS.theta.num/SS.theta.den
            SS.res <- sum(resid(aov.obj)^2)
            SS.res.1 <- SS.res - SS.theta
            T.1df <- ((na * nb - na - nb) * SS.theta)/SS.res.1
            p.value <- 1 - pf(T.1df, 1, na * nb - na - nb)
            list(T.1df = T.1df, p.value = p.value)
      }
 comp.plot_function(aov.obj, data)
{
       vnames <- names(aov.obj$contrasts)
       if(length(vnames) != 2)
             stop("the model must be two-way")
       vara <- data[, vnames[1]]
       varb <- data[, vnames[2]]
       cfs <- coef(aov.obj)
       alpha.A <- aov.obj$contrasts[[vnames[1]]] %*% cfs[
                  aov.obj$assign[[vnames[1]]]]
       alpha.B <- aov.obj$contrasts[[vnames[2]]] %*% cfs[
                  aov.obj$assign[[vnames[2]]]]
       cij <- alpha.B %*% t(alpha.A)
       cij <- c(cij)/cfs[aov.obj$assign$"(Intercept)"]
       na <- length(levels(vara))
       nb <- length(levels(varb))
       r.mat <- matrix(NA, nb, na)
       r.mat[cbind(as.vector(unclass(varb)), as.vector(
                             unclass(vara)))] <- resid(aov.obj)
       plot(cij, as.vector(r.mat))
       ls.fit <- lsfit(as.vector(cij), as.vector(r.mat))
       abline(ls.fit)
       output <- ls.print(ls.fit, print.it = F)
       list(theta.hat = output$coef.table[2, 1], std.error
                      = output$coef.table[2, 2], R.squared =
                        output$summary[2])
}
 tukey.1(pen.aov,pen.df)
 comp.plot(pen.aov,pen.df)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 is.random(pen.df$blend)_T
 is.random(pen.df)

 pen.vc_varcomp(yield~blend,pen.df)
 summary(pen.vc)