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