R/nonadditivity.R

In myaseen208/agricolae: Statistical Procedures for Agricultural Research

#' Nonadditivity model test
#' 
#' The resistance for the transformable nonadditivity, due to J. W. Tukey, is
#' based on the detection of a curvilinear relation between y-est(y) and
#' est(y). A freedom degree for the transformable nonadditivity.
#' 
#' Only two factor: Block and treatment or factor 1 and factor 2.
#' 
#' @param y Answer of the experimental unit
#' @param factor1 Firts treatment applied to each experimental unit
#' @param factor2 Second treatment applied to each experimental unit
#' @param df Degrees of freedom of the experimental error
#' @param MSerror Means square error of the experimental
#' @return P, Q and non-additivity analysis of variance
#' @author Felipe de Mendiburu
#' @references 1. Steel, R.; Torri,J; Dickey, D.(1997) Principles and
#' Procedures of Statistics A Biometrical Approach
#' 
#' 2. George E.P. Box; J. Stuart Hunter and William G. Hunter.  Statistics for
#' experimenters.  Wile Series in probability and statistics
#' @keywords models
#' @importFrom stats residuals coef
#' @export
#' @examples
#' 
#' library(agricolae)
#' data(potato )
#' potato[,1]<-as.factor(potato[,1])
#' model<-lm(cutting ~ date + variety,potato)
#' df<-df.residual(model)
#' MSerror<-deviance(model)/df
#' analysis<-with(potato,nonadditivity(cutting, date, variety, df, MSerror))
#' 
nonadditivity <-
function (y,factor1, factor2, df, MSerror)
{
name.y <- paste(deparse(substitute(y)))
conjunto <- subset(data.frame(y, factor1,factor2), is.na(y) == FALSE)
zz<-tapply(conjunto[,1],conjunto[,c(2,3)],mean)
nn<-tapply(conjunto[,1],conjunto[,c(2,3)],length)
r<-1/mean(1/nn) # Repeticiones
sc<-df*MSerror
x1 <- rownames(zz)
y1 <- colnames(zz)
fila <- length(x1)
col <- length(y1)
total <- fila * col
x <- character(length = total)
y <- character(length = total)
z <- numeric(length = total)
k <- 0
for (i in 1:fila) {
    for (j in 1:col) {
        k <- k + 1
        x[k] <- x1[i]
        y[k] <- y1[j]
        z[k] <- zz[i, j]
    }
}
yy <- data.frame(x, y, z)
yy[,1] <- as.factor(yy[,1])
yy[,2] <- as.factor(yy[,2])
modelo.y<-lm(yy[,3] ~ yy[,1]+ yy[,2])
y.est<-predict(modelo.y)
residual  <-residuals(modelo.y)
q<-y.est^2
modelo.q <- lm(q ~ yy[,1]+ yy[,2])
Nonadditivity  <-residuals(modelo.q)
modelo.e <- lm( residual ~ Nonadditivity )
anva <- anova(modelo.e)
beta <-coef(modelo.e)[2]
P <-  r*anva[1,2] / beta
Q <- P/beta
#one degree free for non-additivity
Sna <-as.numeric(P^2/Q)
if( r == 1) {
sc<-sc-Sna
df<-df-1
}
# Construction of the non-additivity table
anva<-anova(modelo.e)
anva["Residuals",1]<-df
anva["Residuals",2]<-sc
anva["Residuals",3]<-sc/df
anva["Nonadditivity",2:3]<-r*anva["Nonadditivity",2:3]
anva["Nonadditivity",4]<-anva["Nonadditivity",3]/anva["Residuals",3]
anva["Nonadditivity",5]<- 1- pf( anva["Nonadditivity",4],1, df)
cat("\nTukey's test of nonadditivity\n")
cat(name.y,"\n")
cat("\nP :", P)
cat("\nQ :", Q,"\n\n")
print(anva)
output <- list(P=P, Q=Q, ANOVA=anva)
return(output)
}