# R/ea2.R In easyanova: Analysis of Variance and Other Important Complementary Analyses

#### Documented in ea2

ea2<-
function(data, design=1, alpha=0.05, cov=4, list=FALSE, p.adjust=1, plot=2)
{
list=ifelse(list==FALSE,1,2)

pres=function(m){
r=resid(m)
r=scale(r)
t=1:length(r)
g1=function(r){boxplot(r,col="grey80",ylab="Standardized residuals", main="Box plot for residuals")}
g2=function(r){plot(r~t, pch="",ylim=c(-4,4), ylab="Standardized residuals", xlab="Sequence data", 		main="Standardized residuals vs Sequence data", axes=FALSE);axis(2,c(-4,-3.5,-3,-2.5,-2,-1,0,1,2,2.5,3,3.5,4));abline(h=2.5, lty=2);abline(h=-2.5,lty=2);abline(h=3.5, lty=2, col=2);abline(h=-3.5,lty=2, col=2); text(2.5,2.7, "2.5 z-score");text(2.5,-2.7, "-2.5 z-score");text(2.5,3.7, "3.5 z-score");text(2.5,-3.7, "-3.5 z-score");text(t,r,labels=1:length(r))}
a=qqnorm(r,plot.it = FALSE)
a1=a\$x;a2=a\$y;rownames(a2)=NULL; a3=sqrt((a\$y)^2);rownames(a3)=NULL; a4=1:length(r)
d=data.frame(a1,a2,a3,a4)
do=d[order(d[,3], decreasing=TRUE),]
d1=do[1,c(1,2)]
d2=do[2,c(1,2)]
d3=do[3,c(1,2)]
n1=as.character(do[1,4])
n2=as.character(do[2,4])
n3=as.character(do[3,4])
g=list(g1,g2,g3)
g[[plot]](r)
}

sk=function(means, df1, QME, nrep, alpha=0.05){
sk1=function(means, df1, QME, nrep, alpha=alpha) {
means=sort(means,decreasing=TRUE)
n=1:(length(means)-1)
n=as.list(n)
f=function(n){list(means[c(1:n)],means[-c(1:n)])}
g=lapply(n, f)
b1=function(x){(sum(g[[x]][[1]])^2)/length(g[[x]][[1]]) +
(sum(g[[x]][[2]])^2)/length(g[[x]][[2]])-
(sum(c(g[[x]][[1]],g[[x]][[2]]))^2)/length(c(g[[x]][[1]],g[[x]][[2]]))}
p=1:length(g)
values=sapply(p,b1)
minimo=min(values); maximo=max(values)
alfa=(1/(length(means)+df1))*(sum((means-mean(means))^2)+(df1*QME/nrep))
lambda=(pi/(2*(pi-2)))*(maximo/alfa)
vq=qchisq((alpha),lower.tail=FALSE, df=length(means)/(pi-2))
ll=1:length(values); da=data.frame(ll,values); da=da[order(-values),]
ran=da\$ll[1]
r=g[[ran]]; r=as.list(r)
i=ifelse(vq>lambda|length(means)==1, 1,2)
means=list(means)
res=list(means, r)
return(res[[i]])
}
u=sk1(means, df1, QME, nrep, alpha=alpha)
u=lapply(u, sk1, df1=df1, QME=QME, nrep=nrep, alpha=alpha)
sk2=function(u){
v1=function(...){c(u[[1]])};v2=function(...){c(u[[1]],u[[2]])};v3=function(...){c(u[[1]],u[[2]],u[[3]])}
v4=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]])}; v5=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]])}
v6=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]])}
v7=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]])}
v8=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]])}
v9=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]])}
v10=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]])}
v11=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]])}
v12=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]])}
v13=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]])}
v14=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]], u[[14]])}
v15=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]], u[[14]], u[[15]])}
lv=list(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15)
l=length(u)
ti=lv[[l]]
u=ti()
u=lapply(u, sk1, df1=df1, QME=QME, nrep=nrep, alpha=alpha)
return(u)
}
u=sk2(u);u=sk2(u);u=sk2(u);u=sk2(u);u=sk2(u)
u=sk2(u);u=sk2(u);u=sk2(u);u=sk2(u);u=sk2(u)
v1=function(...){c(u[[1]])};v2=function(...){c(u[[1]],u[[2]])};v3=function(...){c(u[[1]],u[[2]],u[[3]])}
v4=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]])}; v5=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]])}
v6=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]])}
v7=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]])}
v8=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]])}
v9=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]])}
v10=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]])}
v11=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]])}
v12=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]])}
v13=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]])}
v14=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]], u[[14]])}
v15=function(...){c(u[[1]],u[[2]],u[[3]],u[[4]],u[[5]],u[[6]],u[[7]],u[[8]],u[[9]],u[[10]],u[[11]], u[[12]], u[[13]], u[[14]], u[[15]])}
lv=list(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15)
l=length(u)
ti=lv[[l]]
u=ti()
rp=u
l2=lapply(rp, length)
l2=unlist(l2)
rp2=rep(letters[1:length(rp)], l2)
return(rp2)
}

cv <- function(x) {
sd = (deviance(x)/df.residual(x))^0.5
mm = mean(fitted(x))
r = 100 * sd/mm
return(round(r, 2))
}

fr=function(m,data){
r=resid(m); names(r)=1:length(r); rp=scale(r)[,1]
s <- shapiro.test(r)
b1<- bartlett.test(r~factor_1, data=data)
b2<- bartlett.test(r~factor_2, data=data)
b3<- bartlett.test(r~treatments, data=data)
cvf=cv(m)
rd=as.data.frame((sort(sqrt(r^2),decreasing=TRUE)))
rl=as.list(rownames(rd))
r1=rl[[1]];r2=rl[[2]];r3=rl[[3]]
d=data.frame(round(s\$"p.value",4),round(b1\$"p.value",4),round(b2\$"p.value",4),round(b3\$"p.value",4), round(cvf,2),as.numeric(r1),as.numeric(r2),as.numeric(r3)); d=t(d)
rownames(d)=c("p.value Shapiro-Wilk test","p.value Bartlett test (factor_1)","p.value Bartlett test (factor_2)","p.value Bartlett test (treatments)","coefficient of variation (%)", "first value most discrepant","second value most discrepant","third value most discrepant")
colnames(d)="values"
l=list("residual analysis"=d,"residuals"=r,"standardized residuals"=rp)

return(l)}

fr2=function(m,data){
data=data.frame(data,nr=c(1:length(data\$response)))
data2<-na.omit(data)
r<-resid(m);  names(r)=1:length(r); rp=scale(r)[,1]
s<-shapiro.test(r)
b1<- bartlett.test(r~plot, data=data2)
b2<- bartlett.test(r~split.plot, data=data2)
b3<- bartlett.test(r~treatments, data=data2)
a1=AIC(m);a1=as.numeric(a1);a1=round(a1,4)
b11=BIC(m);b11=as.numeric(b11);b11=round(b11,4)
rd=as.data.frame((sort(sqrt(r^2),decreasing=TRUE)))
rl=as.list(rownames(rd))
r1=rl[[1]];r2=rl[[2]];r3=rl[[3]]
d=data.frame(round(s\$"p.value",4),round(b1\$"p.value",4),round(b2\$"p.value",4),round(b3\$"p.value",4), a1,b11,as.numeric(r1),as.numeric(r2),as.numeric(r3)); d=t(d)
rownames(d)=c("p.value Shapiro-Wilk test","p.value Bartlett test (plot)","p.value Bartlett test (split.plot)","p.value Bartlett test (plot*split.plot)","AIC","BIC", "first value most discrepant","second value most discrepant","third value most discrepant")
colnames(d)="values"
l=list("residual analysis"=d,"residuals"=r,"standardized residuals"=rp)
return(l)}

fr3=function(m,data){
r=resid(m);  names(r)=1:length(r); rp=scale(r)[,1]
s <- shapiro.test(r)
b1<- bartlett.test(r~factor_1, data=data)
b2<- bartlett.test(r~factor_2, data=data)
b3<- bartlett.test(r~factor_3, data=data)
cvf=cv(m)
rd=as.data.frame((sort(sqrt(r^2),decreasing=TRUE)))
rl=as.list(rownames(rd))
r1=rl[[1]];r2=rl[[2]];r3=rl[[3]]
d=data.frame(round(s\$"p.value",4),round(b1\$"p.value",4),round(b2\$"p.value",4),round(b3\$"p.value",4), round(cvf,2),as.numeric(r1),as.numeric(r2),as.numeric(r3)); d=t(d)
rownames(d)=c("p.value Shapiro-Wilk test","p.value Bartlett test (factor_1)","p.value Bartlett test (factor_2)","p.value Bartlett test (factor_3)","coefficient of variation (%)", "first value most discrepant","second value most discrepant","third value most discrepant")
colnames(d)="values"
l=list("residual analysis"=d,"residuals"=r,"standardized residuals"=rp)
return(l)}

fr4=function(m,data){
data=data.frame(data,nr=c(1:length(data\$response)))
data2<-na.omit(data)
r<-resid(m);names(r)=1:length(r); rp=scale(r)[,1]
s<-shapiro.test(r)
b1<- bartlett.test(r~factor_1, data=data2)
b2<- bartlett.test(r~factor_2, data=data2)
b3<- bartlett.test(r~factor_3, data=data2)
a1=AIC(m);a1=as.numeric(a1);a1=round(a1,4)
b11=BIC(m);b11=as.numeric(b11);b11=round(b11,4)
names(r)=data2\$nr
rd=as.data.frame((sort(sqrt(r^2),decreasing=TRUE)))
rl=as.list(rownames(rd))
r1=rl[[1]];r2=rl[[2]];r3=rl[[3]]
d=data.frame(round(s\$"p.value",4),round(b1\$"p.value",4),round(b2\$"p.value",4),round(b3\$"p.value",4), a1,b11,as.numeric(r1),as.numeric(r2),as.numeric(r3)); d=t(d)
rownames(d)=c("p.value Shapiro-Wilk test","p.value Bartlett test (factor_1)","p.value Bartlett test (factor_2)","p.value Bartlett test (factor_3)","AIC","BIC", "first value most discrepant","second value most discrepant","third value most discrepant")
colnames(d)="values"
l=list("residual analysis"=d,"residuals"=r,"standardized residuals"=rp)
return(l)}

fa2=function(a,m){
ress=c(m\$df.residual,sum(m\$),sum(m\$)/m\$df.residual,NA,NA)
res=a; d=data.frame(res[-1,]); d=data.frame(d[,1],d[,2],d[,2]/d[,1],d[,5],d[,6])
;d=rbind(d,ress);d=round(d,4);d1=d[,5]; d2=ifelse(d1<0.001, "<0.001", d1);
d2=d2[-length(d2)];d2=c(d2,"-"); d=d[,-5];d=data.frame(d,d2);d[is.na(d)] <- "-"
names(d)=c("df", "type III SS", "mean square", "F value", "p>F"); rownames(d)=c(rownames(res[-1,]),"residuals")
return(d)
}

fmm=function(ma,dff){
ma=data.frame(ma,co=ma[,1])
ma=ma[order(ma[,2], decreasing=TRUE),]
j=ma[,1];j=as.character(j)
aux <- combn(j, 2)
w <- apply(aux, 2, paste, collapse = " - ")
jj=ma[,2]
auxj <- combn(jj, 2)
yi=auxj[1,]-auxj[2,]
jjj=ma\$standard.error^2
auxjj <- combn(jjj, 2)
si=sqrt((auxjj[1,]+auxjj[2,])/2)
yx=yi/si; yx=yx^2; yx=sqrt(yx)
nmeans=length(ma[,1])
ft=function(yx, nmeans){1-ptukey(yx,nmeans, dff)}
st=ft(yx,nmeans)
st=round(st,4)
fs=function(ns){s=2:ns;return(s)}
ns=nmeans:2
se=sapply(ns, fs)
ns=unlist(se)
ssnk=ft(yx, ns)
ssnk=round(ssnk,4)
sd=1-ptukey(yx,ns, dff)^(1/(ns-1))
sd=round(sd,4)
yxx=yi/(si*sqrt(2))
vt=1-pt(yxx,dff); vt=vt*2
vt=round(vt,4)
lp=list("none","holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr")
ggs=data.frame(w,round(yi,4),st, ssnk, sd, round(pf,4))
return(ggs)
}

fm=function(ma,dff){
ma=data.frame(ma,co=ma[,1])
ma=ma[order(ma[,2], decreasing=TRUE),]
j=ma[,1];j=as.character(j)
aux <- combn(j, 2)
w <- apply(aux, 2, paste, collapse = " - ")
jj=ma[,2]
auxj <- combn(jj, 2)
yi=auxj[1,]-auxj[2,]
jjj=ma\$sem^2
auxjj <- combn(jjj, 2)
si=sqrt((auxjj[1,]+auxjj[2,])/2)
yx=yi/si; yx=yx^2; yx=sqrt(yx)
nmeans=length(ma[,1])
ftt=function(yx, nmeans){1-ptukey(yx,nmeans, dff)}
st=ftt(yx,nmeans)
st=round(st,4)
fs=function(ns){s=2:ns;return(s)}
ns=nmeans:2
se=sapply(ns, fs)
ns=unlist(se)
ssnk=ftt(yx, ns)
ssnk=round(ssnk,4)
sd=1-ptukey(yx,ns, dff)^(1/(ns-1))
sd=round(sd,4)
yxx=yi/(si*sqrt(2))
vt=1-pt(yxx,dff); vt=vt*2
vt=round(vt,4)
lp=list("none","holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr")
ggs=data.frame(w,round(yi,4),st, ssnk, sd, round(pf,4))
return(ggs)
}
ft=function(test, alpha=0.05){
level=alpha
tes1=test[,3]
tes2=test[,4]
tes3=test[,5]
tes4=test[,6]
names(tes1)=test\$pair
names(tes2)=test\$pair
names(tes3)=test\$pair
names(tes4)=test\$pair
tes1=ifelse(tes1<=level,TRUE,FALSE)
tes2=ifelse(tes2<=level,TRUE,FALSE)
tes3=ifelse(tes3<=level,TRUE,FALSE)
tes4=ifelse(tes4<=level,TRUE,FALSE)
x1=tes1;x2=tes2;x3=tes3;x4=tes4
inab <- function(x, Letters=c(letters, LETTERS), separator=".", decreasing = decreasing){
obj_x <- deparse(substitute(x))
namx <- names(x)
namx <- gsub(" ", "", names(x))
if(length(namx) != length(x))
stop("Names required for ", obj_x)
split_names <- strsplit(namx, "-")
stopifnot( sapply(split_names, length) == 2 )
comps <- t(as.matrix(as.data.frame(split_names)))
rownames(comps) <- names(x)
lvls <- unique(as.vector(comps))
n <- length(lvls)
lmat <- array(TRUE, dim=c(n,1), dimnames=list(lvls, NULL) )
if( sum(x) == 0 ){
ltrs <- rep(get_letters(1, Letters=Letters, separator=separator), length(lvls) )
names(ltrs) <- lvls
colnames(lmat) <- ltrs[1]
msl <- ltrs
ret <- list(Letters=ltrs, monospacedLetters=msl, LetterMatrix=lmat)
return(ret)
}
else{
signifs <- comps[x,,drop=FALSE]
absorb <- function(m){
for(j in 1:(ncol(m)-1)){
for(k in (j+1):ncol(m)){
if( all(m[which(m[,k]),k] & m[which(m[,k]),j]) ){
m <- m[,-k, drop=FALSE]
return(absorb(m))
}
else if( all(m[which(m[,j]),k] & m[which(m[,j]),j]) ){
m <- m[,-j, drop=FALSE]
return(absorb(m))
}
}
}
return(m)
}
for( i in 1:nrow(signifs) ){
tmpcomp <- signifs[i,]
wassert <- which(lmat[tmpcomp[1],] & lmat[tmpcomp[2],])
if(any(wassert)){
tmpcols <- lmat[,wassert,drop=FALSE]
tmpcols[tmpcomp[2],] <- FALSE
lmat[tmpcomp[1],wassert] <- FALSE
lmat <- cbind(lmat, tmpcols)
colnames(lmat) <- get_letters( ncol(lmat), Letters=Letters,
separator=separator)
if(ncol(lmat) > 1){
lmat <- absorb(lmat)
colnames(lmat) <- get_letters( ncol(lmat),  Letters=Letters,
separator=separator )
}
}
}
}
lmat <- lmat[,order(apply(lmat, 2, sum))]
lmat <- sweepLetters(lmat)
lmat <- lmat[,names(sort(apply(lmat,2, function(x) return(min(which(x))))))]
colnames(lmat) <- get_letters( ncol(lmat),  Letters=Letters,
separator=separator)
lmat <- lmat[,order(apply(lmat, 2, sum))]
lmat <- sweepLetters(lmat)
lmat <- lmat[,names(sort(apply(lmat,2, function(x) return(min(which(x)))),
decreasing = decreasing))]
colnames(lmat) <- get_letters( ncol(lmat),  Letters=Letters,
separator=separator)
ltrs <- apply(lmat,1,function(x) return(paste(names(x)[which(x)], sep="", collapse="") ) )
msl <- matrix(ncol=ncol(lmat), nrow=nrow(lmat))
for( i in 1:nrow(lmat) ){
msl[i,which(lmat[i,])] <- colnames(lmat)[which(lmat[i,])]
absent <- which(!lmat[i,])
if( length(absent) < 2 ){
if( length(absent) == 0 )
next
else{
msl[i,absent] <- paste( rep(" ", nchar(colnames(lmat)[absent])), collapse="" )
}
}
else{
msl[i,absent] <- unlist( lapply( sapply( nchar(colnames(lmat)[absent]),
function(x) return(rep( " ",x)) ),
paste, collapse="") )
}
}
msl <- apply(msl, 1, paste, collapse="")
names(msl) <- rownames(lmat)
ret <- list(Letters=ltrs)
return(ret)
}

sweepLetters <- function(mat, start.col=1, Letters=c(letters, LETTERS), separator="."){
stopifnot( all(start.col %in% 1:ncol(mat)) )
locked <- matrix(rep(0,ncol(mat)*nrow(mat)), ncol=ncol(mat))
cols <- 1:ncol(mat)
cols <- cols[c( start.col, cols[-start.col] )]
if( any(is.na(cols) ) )
cols <- cols[-which(is.na(cols))]
for( i in cols){
tmp <- matrix(rep(0,ncol(mat)*nrow(mat)), ncol=ncol(mat))
tmp[which(mat[,i]),] <- mat[which(mat[,i]),]
one <- which(tmp[,i]==1)
if( all(apply(tmp[,-i,drop=FALSE], 1, function(x) return( any(x==1) ))) ){
}
for( j in one ){
if( locked[j,i] == 1 ){
next
}
chck <- 0
lck <- list()
for( k in one ){
if( j==k ){
next
}
else{
rows <- tmp[c(j,k),]
dbl <- rows[1,] & rows[2,]
hit <- which(dbl)
hit <- hit[-which(hit==i)]
dbl <- rows[1,-i,drop=FALSE] & rows[2,-i,drop=FALSE]
if( any(dbl) ){
chck <- chck + 1
lck[[chck]] <- list(c(j,hit[length(hit)]), c(k,hit[length(hit)]))
}
}
}
if( (chck == (length(one)-1)) && chck != 0 ){
for( k in 1:length(lck) ){
locked[ lck[[k]][[1]][1], lck[[k]][[1]][2] ] <- 1
locked[ lck[[k]][[2]][1], lck[[k]][[2]][2] ] <- 1
}
mat[j,i] <- FALSE
}
}
if(all(mat[,i]==FALSE)){
mat <- mat[,-i,drop=FALSE]
colnames(mat) <- get_letters( ncol(mat), Letters=Letters, separator=separator)
return(sweepLetters(mat, Letters=Letters, separator=separator))
}
}
onlyF <- apply(mat, 2, function(x) return(all(!x)))
if( any(onlyF) ){
mat <- mat[,-which(onlyF),drop=FALSE]
colnames(mat) <- get_letters( ncol(mat), Letters=Letters, separator=separator)
}
return( mat )
}

get_letters <- function( n, Letters=c(letters, LETTERS), separator="." ){
n.complete <- floor(n / length(Letters))
n.partial <- n %% length(Letters)
lett <- character()
separ=""
if( n.complete > 0 ){
for( i in 1:n.complete ){
lett <- c(lett, paste(separ, Letters, sep="") )
separ <- paste( separ, separator, sep="" )
}
}
if(n.partial > 0 )
lett <- c(lett, paste(separ, Letters[1:n.partial], sep="") )
return(lett)
}
decreasing=FALSE;jjj1=inab(x1, decreasing = decreasing,); jjj1=jjj1[[1]]
jjj2=inab(x2, decreasing = decreasing,); jjj2=jjj2[[1]]
jjj3=inab(x3, decreasing = decreasing,); jjj3=jjj3[[1]]
jjj4=inab(x4, decreasing = decreasing,); jjj4=jjj4[[1]]
return(hgy)
}

fsd1=function(data)
{
s=split(data,data\$factor_1)
f=function(x){sd(s[[x]]\$response, na.rm=TRUE)}
sds=sapply(1:length(s), f); sds=round(sds,4)
return(sds)
}

fsd2=function(data)
{
s=split(data,data\$factor_2)
f=function(x){sd(s[[x]]\$response, na.rm=TRUE)}
sds=sapply(1:length(s), f); sds=round(sds,4)
return(sds)
}

fsd3=function(data)
{
s=split(data,data\$treatment)
f=function(x){sd(s[[x]]\$response, na.rm=TRUE)}
sds=sapply(1:length(s), f); sds=round(sds,4)
return(sds)
}

f1<-function(data,cov){
names(data)=c("factor_1","factor_2","response")
data<-data.frame(factor_1=factor(data\$factor_1), factor_2=factor(data\$factor_2), response=data\$response)
m<-aov(response~factor_1*factor_2,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum))
m1<-aov(response~-1+factor_1+factor_2+factor_1*factor_2,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum))
m2<-aov(response~-1+ factor_2+factor_1+factor_2*factor_1,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum))
a<-anova(m)
a2 <- drop1(m,.~.,test="F")
a3 <- a2
a3<-fa2(a3,m)
treatments=interaction(data\$factor_1,data\$factor_2)
data<-data.frame(data,treatments)
m3<-aov(response~-1+treatments,data=data, contrasts=list(treatments=contr.sum))
data2<-na.omit(data)
res=fr(m,data2)

sem<-round(sderrs(m1)[1:nlevels(data\$factor_1)],4)
sd=fsd1(data)
factor_1<-levels(data\$factor_1)
rownames(maf1)=NULL

sem<-round(sderrs(m2) [1:nlevels(data\$factor_2)],4)
factor_2<-levels(data\$factor_2)
sd=fsd2(data)
rownames(maf2)=NULL

sem<-round(sderrs(m3),4)
treatment<-levels(data\$t)
sd=fsd3(data)
rownames(mat)=NULL

dff=df.residual(m)
test1=fm(maf1,dff)
test2=fm(maf2,dff)
nrep1=length(,1])/length(factor_1)
nrep2=length(,1])/length(factor_2)
nrep3=length(,1])/length(treatment)
QME=deviance(m)/dff
scott_knott=sk(means1, dff, QME, nrep1, alpha)
groups1=ft(test1, alpha); maf1=maf1[order(maf1[,2], decreasing=TRUE),]
mf1=data.frame(maf1,groups1, scott_knott)
rownames(mf1) = NULL
groups2=ft(test2, alpha); maf2=maf2[order(maf2[,2], decreasing=TRUE),]
scott_knott=sk(means2, dff, QME, nrep2, alpha)
mf2=data.frame(maf2,groups2, scott_knott)
rownames(mf2) = NULL
c1=rep(1:nlevels(data\$factor_2), each=nlevels(data\$factor_1))
cs1=split(mat, c1)
test3=lapply(cs1, fm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_2));n11=data.frame(n1,levels(data\$factor_2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2))
cs2=split(mat, c2)
test4=lapply(cs2, fm, dff=dff); n2=rep("factor_2 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf)
i=1:length(mft1)
ft1=function(i){a=mft1[[i]][,2];names(a)=mft1[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups3[x],scott_knott[x])}
xx1=1:length(groups3)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups3)
dd<-function(x){l=data.frame(mft1[[x]],ftr[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_2)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf)
i=1:length(mft2)
ft1=function(i){a=mft2[[i]][,2];names(a)=mft2[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups4[x],scott_knott[x])}
xx1=1:length(groups4)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups4)
ddd<-function(x){l=data.frame(mft2[[x]],ftr[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
l<-list(a3,mf1, test1, mf2,test2,mft1,test3, mft2, test4, res)
names(l)= list("Analysis of variance", "Adjusted means (factor 1)", "Multiple comparison test (factor 1)","Adjusted means (factor 2)", "Multiple comparison test (factor 2)", "Adjusted means (factor 1 in levels of factor 2)", "Multiple comparison test (factor 1 in levels of factor 2)", "Adjusted means (factor 2 in levels of factor 1)", "Multiple comparison test (factor 2 in levels of factor 1)","Residual analysis")
pres(m)
return(l)}

f2<-function(data, cov){
names(data)=c("factor_1","factor_2","blocks","response")
data<-data.frame(factor_1=factor(data\$factor_1), factor_2=factor(data\$factor_2), blocks=factor(data\$blocks), response=data\$response)
m<-aov(response~factor_1*factor_2+blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, blocks=contr.sum))
m1<-aov(response~-1+factor_1+factor_2+factor_1*factor_2+blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, blocks=contr.sum))
m2<-aov(response~-1+ factor_2+factor_1+factor_2*factor_1+blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, blocks=contr.sum))
a<-anova(m)
a2 <- drop1(m,.~.,test="F")
a3 <- a2
a3<-fa2(a3,m)
treatments=interaction(data\$factor_1,data\$factor_2)
data<-data.frame(data,treatments)
m3<-aov(response~-1+treatments+blocks,data=data, contrasts=list(treatments=contr.sum, blocks=contr.sum))
data2<-na.omit(data)
res=fr(m,data2)

sem<-round(sderrs(m1)[1:nlevels(data\$factor_1)],4)
sd=fsd1(data)
factor_1<-levels(data\$factor_1)
rownames(maf1)=NULL

sem<-round(sderrs(m2)[1:nlevels(data\$factor_2)],4)
sd=fsd2(data)
factor_2<-levels(data\$factor_2)
rownames(maf2)=NULL

treatment<-levels(data\$treatments)
sem<-round(sderrs(m3)[1:nlevels(data\$treatments)],4)
sd=fsd3(data)
rownames(mat)=NULL

dff=df.residual(m)
test1=fm(maf1,dff)
test2=fm(maf2,dff)
nrep1=length(,1])/length(factor_1)
nrep2=length(,1])/length(factor_2)
nrep3=length(,1])/length(treatment)
QME=deviance(m)/dff
scott_knott=sk(means1, dff, QME, nrep1, alpha)
groups1=ft(test1, alpha); maf1=maf1[order(maf1[,2], decreasing=TRUE),]
mf1=data.frame(maf1,groups1, scott_knott)
rownames(mf1) = NULL
groups2=ft(test2, alpha); maf2=maf2[order(maf2[,2], decreasing=TRUE),]
scott_knott=sk(means2, dff, QME, nrep2, alpha)
mf2=data.frame(maf2,groups2, scott_knott)
rownames(mf2) = NULL
c1=rep(1:nlevels(data\$factor_2), each=nlevels(data\$factor_1))
cs1=split(mat, c1)
test3=lapply(cs1, fm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_2));n11=data.frame(n1,levels(data\$factor_2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2))
cs2=split(mat, c2)
test4=lapply(cs2, fm, dff=dff); n2=rep("factor_2 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf)
i=1:length(mft1)
ft1=function(i){a=mft1[[i]][,2];names(a)=mft1[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups3[x],scott_knott[x])}
xx1=1:length(groups3)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups3)
dd<-function(x){l=data.frame(mft1[[x]],ftr[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_2)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf)
i=1:length(mft2)
ft1=function(i){a=mft2[[i]][,2];names(a)=mft2[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups4[x],scott_knott[x])}
xx1=1:length(groups4)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups4)
ddd<-function(x){l=data.frame(mft2[[x]],ftr[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
l<-list(a3,mf1, test1, mf2,test2,mft1,test3, mft2, test4, res)
names(l)= list("Analysis of variance", "Adjusted means (factor 1)", "Multiple comparison test (factor 1)","Adjusted means (factor 2)", "Multiple comparison test (factor 2)", "Adjusted means (factor 1 in levels of factor 2)", "Multiple comparison test (factor 1 in levels of factor 2)", "Adjusted means (factor 2 in levels of factor 1)", "Multiple comparison test (factor 2 in levels of factor 1)","Residual analysis")
pres(m)
return(l)
}

f3<-function(data, cov){
names(data)=c("factor_1","factor_2","rows","cols","response")
data<-data.frame(factor_1=factor(data\$factor_1), factor_2=factor(data\$factor_2), rows=factor(data\$rows), columns=factor(data\$cols), response=data\$response)
m<-aov(response~factor_1*factor_2+ rows+columns,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, rows=contr.sum, columns=contr.sum))
m1<-aov(response~-1+factor_1+factor_2+factor_1*factor_2+ rows+columns,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, rows=contr.sum, columns=contr.sum))
m2<-aov(response~-1+ factor_2+factor_1+factor_2*factor_1+ rows+columns,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, rows=contr.sum, columns=contr.sum))
a<-anova(m)
a2 <- drop1(m,.~.,test="F")
a3 <- a2
a3<-fa2(a3,m)
treatments=interaction(data\$factor_1,data\$factor_2)
data<-data.frame(data,treatments)
m3<-aov(response~-1+treatments +rows+columns,data=data, contrasts=list(treatments=contr.sum, rows=contr.sum, columns=contr.sum))
data2<-na.omit(data)
res=fr(m,data2)

sem<-round(sderrs(m1)[1:nlevels(data\$factor_1)],4)
sd=fsd1(data)
factor_1<-levels(data\$factor_1)
rownames(maf1)=NULL

sem<-round(sderrs(m2)[1:nlevels(data\$factor_2)],4)
sd=fsd2(data)
factor_2<-levels(data\$factor_2)
rownames(maf2)=NULL

treatment<-levels(data\$treatments)
sem<-round(sderrs(m3)[1:nlevels(data\$treatments)],4)
sd=fsd3(data)
rownames(mat)=NULL

dff=df.residual(m)
nrep1=length(,1])/length(factor_1)
nrep2=length(,1])/length(factor_2)
nrep3=length(,1])/length(treatment)
QME=deviance(m)/dff
scott_knott=sk(means1, dff, QME, nrep1, alpha)
test1=fm(maf1,dff)
test2=fm(maf2,dff)
groups1=ft(test1, alpha); maf1=maf1[order(maf1[,2], decreasing=TRUE),]
mf1=data.frame(maf1,groups1, scott_knott)
rownames(mf1) = NULL
groups2=ft(test2, alpha); maf2=maf2[order(maf2[,2], decreasing=TRUE),]
scott_knott=sk(means2, dff, QME, nrep1, alpha)
mf2=data.frame(maf2,groups2, scott_knott)
rownames(mf2) = NULL
c1=rep(1:nlevels(data\$factor_2), each=nlevels(data\$factor_1))
cs1=split(mat, c1)
test3=lapply(cs1, fm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_2));n11=data.frame(n1,levels(data\$factor_2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2))
cs2=split(mat, c2)
test4=lapply(cs2, fm, dff=dff); n2=rep("factor_2 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf)
i=1:length(mft1)
ft1=function(i){a=mft1[[i]][,2];names(a)=mft1[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups3[x],scott_knott[x])}
xx1=1:length(groups3)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups3)
dd<-function(x){l=data.frame(mft1[[x]],ftr[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_2)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf)
i=1:length(mft2)
ft1=function(i){a=mft2[[i]][,2];names(a)=mft2[[i]][,1];return(a)}
lap1=lapply(i, ft1)
scott_knott=lapply(lap1, sk, df1=dff, QME=QME, nrep=nrep3, alpha=alpha)
names(scott_knott)=rep("scott_knott",length(scott_knott))
ftt1=function(x){a=data.frame(groups4[x],scott_knott[x])}
xx1=1:length(groups4)
ftr=lapply(xx1, ftt1)
names(ftr)=names(groups4)
ddd<-function(x){l=data.frame(mft2[[x]],ftr[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
l<-list(a3,mf1, test1, mf2,test2,mft1,test3, mft2, test4, res)
names(l)= list("Analysis of variance", "Adjusted means (factor 1)", "Multiple comparison test (factor 1)","Adjusted means (factor 2)", "Multiple comparison test (factor 2)", "Adjusted means (factor 1 in levels of factor 2)", "Multiple comparison test (factor 1 in levels of factor 2)", "Adjusted means (factor 2 in levels of factor 1)", "Multiple comparison test (factor 2 in levels of factor 1)","Residual analysis")
pres(m)
return(l)
}

f4<-function(data, cov){
gg<-cov
names(data)<-c("plot","rep","split.plot","response")
data<-data.frame(plot=factor(data\$plot), rep=factor(data\$rep), split.plot=factor(data\$split.plot), response=data\$response)
subject<-interaction(data\$plot,data\$rep)
treatments<-interaction(data\$plot,data\$split.plot)
data<-data.frame(data, subject,treatments)
UN<-corSymm(form=~1|subject)
AR<-corAR1(form=~1|subject)
ARH<-corAR1(form=~1|subject)
CS<- corCompSymm (form=~1|subject)
CAR<-corCAR1(form=~split.plot|subject)
o1<-list(AR,ARH, CAR, CS, UN)
cor<-o1[[gg]]
UN1<-varIdent(form=~1|split.plot)
AR1<-NULL
ARH1<- varIdent(form=~1|split.plot)
CS1<-NULL
CAR1<-NULL
o2<-list(AR1,ARH1, CAR1, CS1, UN1)
var<-o2[[gg]]
m<-lme(response~plot*split.plot, random=~1|subject, correlation=cor, weights=var, data<-data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum), control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
b<-anova(m, type="marginal")[-1,]
m1<-lme(response~-1+ treatments, random=~1|subject, data=data, na.action=na.omit, contrasts=list(treatments=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
am<-fixef(m1)[1:nlevels(data\$treatments)]
c<-data.frame(levels(treatments), round(am,4),
round(sderrs(m1)[1:nlevels(data\$treatments)],4))
rownames(c)<-NULL
df1<- anova(m) \$"denDF"[2]
df2<- anova(m) \$"denDF"[1]
mm1<-lme(response~-1+plot*split.plot, random=~1|subject, data=data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000,  niterEM=2000, opt="optim"))
mm2<-lme(response~-1+split.plot*plot, random=~1|subject, data=data, na.action=na.omit, contrasts=list(split.plot=contr.sum, plot=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
a11<-data.frame(levels(data\$plot), round(fixef(mm1)[1:nlevels(data\$plot)],4), round(sderrs(mm1)[1:nlevels(data\$plot)],4))
a22<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(sderrs(mm2)[1:nlevels(data\$split.plot)],4))

#################################################
# MS
msa=sqrt((m\$sigma^2)/(nlevels(,1])*(nrow(data)/nlevels(data\$treatments))))
msa=rep(msa,nlevels(,3]))

a222<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(msa,4))

a22<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(sderrs(mm2)[1:nlevels(data\$split.plot)],4))

#################################################
mal=lm(response~plot+subject, data=data)
ano=anova(mal)
vara=anova(mal)[[2,3]]
varb=m\$sigma^2
J=nlevels(data\$split.plot)-1
df1a=((vara+(J*varb))^2)/( ((vara^2)/df1)+(((J*varb)^2)/df2))

rownames(a11)<-NULL
rownames(a22)<-NULL
rownames(a222)<-NULL
res=fr2(m,data)
test1=fmm(a11,df1)
test2=fmm(a222,df2)
groups1=ft(test1, alpha); a11=a11[order(a11[,2], decreasing=TRUE),]
mf1=data.frame(a11,groups1);rownames(mf1) = NULL
groups2=ft(test2, alpha); a22=a22[order(a22[,2], decreasing=TRUE),]
mf2=data.frame(a22,groups2);rownames(mf2) = NULL
c1=rep(1:nlevels(data\$split.plot), each=nlevels(data\$plot))
cs1=split(c, c1)
test3=lapply(cs1, fmm, dff=df1a); n1=rep("plot in", nlevels(data\$split.plot));n11=data.frame(n1,levels(data\$split.plot))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$plot), nlevels(data\$split.plot))
cs2=split(c, c2)
###############
msa2=sqrt((m\$sigma^2)/(nrow(data)/nlevels(data\$treatments)))
msa2=rep(msa2,nlevels(data\$treatments))
c222<-data.frame(levels(treatments), round(am,4), round(msa2,4))
rownames(c222)<-NULL
cs222=split(c222, c2)

#######################
test4=lapply(cs222, fmm, dff=df2); n2=rep("split.plot in", nlevels(data\$plot));n22=data.frame(n2,levels(data\$plot))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft1[[x]],groups3[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$split.plot)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft2[[x]],groups4[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$plot)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
cva=round(sqrt(vara)*100/mean(data\$response, na.rm = TRUE),4); vara=round(vara,4)
cvb=round(sqrt(varb)*100/mean(data\$response, na.rm = TRUE),4);varb=round(varb,4)

ress=rbind(res[[1]],"Mean Square of Error a"=vara, "Mean Square of Error b"=varb, "Coefficient of Variation a"=cva,"Coefficient of Variation b"=cvb)
res=list("values"=ress,"residuals"=res[[2]],"standardized residuals"=res[[3]])

l=list(b,mf1,test1,mf2,test2,mft1,test3,mft2,test4, res)
names(l)=c("Marginal anova (Type III Sum of Squares)","Adjusted means (plot)", "Multiple comparison test (plot)","Adjusted means (split.plot)", "Multiple comparison test (split.plot)", "Adjusted means (plot in levels of split.plot)", "Multiple comparison test (plot in levels of split.plot)", "Adjusted means (split.plot in levels of plot)", "Multiple comparison test (split.plot in levels of plot)","Residual analysis")
pres(m)
return(l)
}

f5<-function(data, cov){
gg<-cov
names(data)<-c("plot","block","split.plot","response")
data<-data.frame(plot=factor(data\$plot), block=factor(data\$block), split.plot=factor(data\$split.plot), response=data\$response)
subject<-interaction(data\$plot,data\$block)
treatments<-interaction(data\$plot,data\$split.plot)
data<-data.frame(data, subject,treatments)
UN<-corSymm(form=~1|subject)
AR<-corAR1(form=~1|subject)
ARH<-corAR1(form=~1|subject)
CS<- corCompSymm (form=~1|subject)
CAR<-corCAR1(form=~split.plot|subject)
o1<-list(AR,ARH, CAR, CS, UN)
cor<-o1[[gg]]
UN1<-varIdent(form=~1|split.plot)
AR1<-NULL
ARH1<- varIdent(form=~1|split.plot)
CS1<-NULL
CAR1<-NULL
o2<-list(AR1,ARH1, CAR1, CS1, UN1)
var<-o2[[gg]]
m<-lme(response~plot*split.plot+block, random=~1|subject, correlation=cor, weights=var, data<-data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum, block=contr.sum), control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
b<-anova(m, type="marginal")[-1,]
m1<-lme(response~-1+ treatments+block, random=~1|subject, data=data, na.action=na.omit, contrasts=list(treatments=contr.sum, block=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
am<-fixef(m1)[1:nlevels(data\$ treatments)]
c<-data.frame(levels(treatments), round(am,4), round(sderrs(m1)[1:nlevels(data\$ treatments)],4))
rownames(c)<-NULL
df1<- anova(m) \$"denDF"[2]
df2<- anova(m) \$"denDF"[1]
mm1<-lme(response~-1+plot*split.plot+block, random=~1|subject, data=data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum, block=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
mm2<-lme(response~-1+split.plot*plot+block, random=~1|subject, data=data, na.action=na.omit, contrasts=list(split.plot=contr.sum, plot=contr.sum, block=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
a11<-data.frame(levels(data\$plot), round(fixef(mm1)[1:nlevels(data\$plot)],4), round(sderrs(mm1)[1:nlevels(data\$plot)],4))
#################################################
# MS
msa=sqrt((m\$sigma^2)/(nlevels(,1])*(nrow(data)/nlevels(data\$treatments))))
msa=rep(msa,nlevels(,3]))

a222<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(msa,4))

a22<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(sderrs(mm2)[1:nlevels(data\$split.plot)],4))

#################################################
mal=lm(response~block*plot, data=data)
ano=anova(mal)
vara=anova(mal)[[3,3]]
varb=m\$sigma^2

J=nlevels(data\$split.plot)-1
df1a=((vara+(J*varb))^2)/( ((vara^2)/df1)+(((J*varb)^2)/df2))

rownames(a11)<-NULL
rownames(a22)<-NULL
rownames(a222)<-NULL
res=fr2(m,data)
test1=fmm(a11,df1)
test2=fmm(a222,df2)
groups1=ft(test1, alpha); a11=a11[order(a11[,2], decreasing=TRUE),]
mf1=data.frame(a11,groups1);rownames(mf1) = NULL
groups2=ft(test2, alpha); a22=a22[order(a22[,2], decreasing=TRUE),]
mf2=data.frame(a22,groups2);rownames(mf2) = NULL
c1=rep(1:nlevels(data\$split.plot), each=nlevels(data\$plot))
cs1=split(c, c1)
test3=lapply(cs1, fmm, dff=df1a); n1=rep("plot in", nlevels(data\$split.plot));n11=data.frame(n1,levels(data\$split.plot))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$plot), nlevels(data\$split.plot))
cs2=split(c, c2)
###############
msa2=sqrt((m\$sigma^2)/(nrow(data)/nlevels(data\$treatments)))
msa2=rep(msa2,nlevels(data\$treatments))
c222<-data.frame(levels(treatments), round(am,4), round(msa2,4))
rownames(c222)<-NULL
cs222=split(c222, c2)

#######################
test4=lapply(cs222, fmm, dff=df2); n2=rep("split.plot in", nlevels(data\$plot));n22=data.frame(n2,levels(data\$plot))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft1[[x]],groups3[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$split.plot)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft2[[x]],groups4[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$plot)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
cva=round(sqrt(vara)*100/mean(data\$response, na.rm = TRUE),4); vara=round(vara,4)
cvb=round(sqrt(varb)*100/mean(data\$response, na.rm = TRUE),4);varb=round(varb,4)

ress=rbind(res[[1]],"Mean Square of Error a"=vara, "Mean Square of Error b"=varb, "Coefficient of Variation a"=cva,"Coefficient of Variation b"=cvb)
res=list("values"=ress,"residuals"=res[[2]],"standardized residuals"=res[[3]])

l=list(b,mf1,test1,mf2,test2,mft1,test3,mft2,test4, res)
names(l)=c("Marginal anova (Type III Sum of Squares)","Adjusted means (plot)", "Multiple comparison test (plot)","Adjusted means (split.plot)", "Multiple comparison test (split.plot)", "Adjusted means (plot in levels of split.plot)", "Multiple comparison test (plot in levels of split.plot)", "Adjusted means (split.plot in levels of plot)", "Multiple comparison test (split.plot in levels of plot)","Residual analysis")
pres(m)
return(l)
}

f6<-function(data, cov){
gg<-cov
names(data)<-c("plot","row", "column", "split.plot","response")
data<-data.frame(plot=factor(data\$plot), row=factor(data\$row), column=factor(data\$column), split.plot=factor(data\$split.plot), response=data\$response)
subject<-interaction(data\$plot,data\$column)
treatments<-interaction(data\$plot,data\$split.plot)
data<-data.frame(data, subject,treatments)
UN<-corSymm(form=~1|subject)
AR<-corAR1(form=~1|subject)
ARH<-corAR1(form=~1|subject)
CS<- corCompSymm (form=~1|subject)
CAR<-corCAR1(form=~split.plot|subject)
o1<-list(AR,ARH, CAR, CS, UN)
cor<-o1[[gg]]
UN1<-varIdent(form=~1|split.plot)
AR1<-NULL
ARH1<- varIdent(form=~1|split.plot)
CS1<-NULL
CAR1<-NULL
o2<-list(AR1,ARH1, CAR1, CS1, UN1)
var<-o2[[gg]]
m<-lme(response~plot*split.plot+row+column, random=~1|subject, correlation=cor, weights=var, data<-data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum, column=contr.sum, row=contr.sum), control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
b<-anova(m, type="marginal")[-1,]
m1<-lme(response~-1+ treatments+row+column, random=~1|subject, data=data, na.action=na.omit, contrasts=list(treatments=contr.sum, column=contr.sum, row=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000,  opt="optim"))
am<-fixef(m1)[1:nlevels(data\$treatments)]
c<-data.frame(levels(treatments), round(am,4), round(sderrs(m1)[1:nlevels(data\$ treatments)],4))
rownames(c)<-NULL
df1<- anova(m) \$"denDF"[2]
df2<- anova(m) \$"denDF"[1]
mm1<-lme(response~-1+plot*split.plot+row+column, random=~1|subject, data=data, na.action=na.omit, contrasts=list(plot=contr.sum, split.plot=contr.sum, column=contr.sum, row=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
mm2<-lme(response~-1+split.plot*plot+row+column, random=~1|subject, data=data, na.action=na.omit, contrasts=list(split.plot=contr.sum, plot=contr.sum, column=contr.sum, row=contr.sum), correlation=cor, weights=var, control=lmeControl(maxIter =6000, msMaxIter=6000, niterEM=2000, opt="optim"))
a11<-data.frame(levels(data\$plot),
round( fixef(mm1)[1:nlevels(data\$plot)], 4),
round(sderrs(mm1)[1:nlevels(data\$plot)], 4))
a22<-data.frame(levels(data\$split.plot),
round( fixef(mm2)[1:nlevels(data\$split.plot)],4),
round(sderrs(mm2)[1:nlevels(data\$split.plot)],4))

#################################################
# MS
msa=sqrt((m\$sigma^2)/(nlevels(,1])*(nrow(data)/nlevels(data\$treatments))))
msa=rep(msa,nlevels(,4]))

a222<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(msa,4))

a22<-data.frame(levels(data\$split.plot), round(fixef(mm2)[1:nlevels(data\$split.plot)],4), round(sderrs(mm2)[1:nlevels(data\$split.plot)],4))

#################################################
mal=lm(response~row+column+plot+subject, data=data)
ano=anova(mal)
vara=anova(mal)[[4,3]]
varb=m\$sigma^2
J=nlevels(data\$split.plot)-1
df1a=((vara+(J*varb))^2)/( ((vara^2)/df1)+(((J*varb)^2)/df2))

rownames(a11)<-NULL
rownames(a22)<-NULL
rownames(a222)<-NULL
res=fr2(m,data)
test1=fmm(a11,df1)
test2=fmm(a222,df2)
groups1=ft(test1, alpha); a11=a11[order(a11[,2], decreasing=TRUE),]
mf1=data.frame(a11,groups1);rownames(mf1) = NULL
groups2=ft(test2, alpha); a22=a22[order(a22[,2], decreasing=TRUE),]
mf2=data.frame(a22,groups2);rownames(mf2) = NULL
c1=rep(1:nlevels(data\$split.plot), each=nlevels(data\$plot))
cs1=split(c, c1)
test3=lapply(cs1, fmm, dff=df1a); n1=rep("plot in", nlevels(data\$split.plot));n11=data.frame(n1,levels(data\$split.plot))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test3)=n11
c2=rep(1:nlevels(data\$plot), nlevels(data\$split.plot))
cs2=split(c, c2)
###############
msa2=sqrt((m\$sigma^2)/(nrow(data)/nlevels(data\$treatments)))
msa2=rep(msa2,nlevels(data\$treatments))
c222<-data.frame(levels(treatments), round(am,4), round(msa2,4))
rownames(c222)<-NULL
cs222=split(c222, c2)

#######################
test4=lapply(cs222, fmm, dff=df2); n2=rep("split.plot in", nlevels(data\$plot));n22=data.frame(n2,levels(data\$plot))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test4)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups3=lapply(test3, ft, alpha)
mft1=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft1[[x]],groups3[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$split.plot)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups4=lapply(test4, ft, alpha)
mft2=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft2[[x]],groups4[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$plot)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
cva=round(sqrt(vara)*100/mean(data\$response, na.rm = TRUE),4); vara=round(vara,4)
cvb=round(sqrt(varb)*100/mean(data\$response, na.rm = TRUE),4);varb=round(varb,4)

ress=rbind(res[[1]],"Mean Square of Error a"=vara, "Mean Square of Error b"=varb, "Coefficient of Variation a"=cva,"Coefficient of Variation b"=cvb)
res=list("values"=ress,"residuals"=res[[2]],"standardized residuals"=res[[3]])

l=list(b,mf1,test1,mf2,test2,mft1,test3,mft2,test4, res)
names(l)=c("Marginal anova (Type III Sum of Squares)","Adjusted means (plot)", "Multiple comparison test (plot)","Adjusted means (split.plot)", "Multiple comparison test (split.plot)", "Adjusted means (plot in levels of split.plot)", "Multiple comparison test (plot in levels of split.plot)", "Adjusted means (split.plot in levels of plot)", "Multiple comparison test (split.plot in levels of plot)","Residual analysis")
pres(m)
return(l)
}

f7<-function(data, cov){
names(data)=c("factor_1","factor_2", "factor_3", "response")
data<-data.frame(factor_1=factor(data\$factor_1), factor_2=factor(data\$factor_2), factor_3=factor(data\$factor_3), response=data\$response)
m<-aov(response~factor_1*factor_2*factor_3,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum))
m1<-aov(response~-1+ factor_1*factor_2*factor_3,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum))
m2<-aov(response~-1+ factor_2*factor_1*factor_3,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum))
m3<-aov(response~-1+ factor_3*factor_1*factor_2,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum))
a2 <- drop1(m,.~.,test="F")
a3 <- a2
a3<-fa2(a3,m)
data2<-na.omit(data)
res=fr3(m,data2)
treatments_f1f2=interaction(data\$factor_1,data\$factor_2)
treatments_f1f3=interaction(data\$factor_1,data\$factor_3)
treatments_f2f3=interaction(data\$factor_2,data\$factor_3)
treatments_f1f2f3=interaction(data\$factor_1,data\$factor_2, data\$factor_3)
treatments_f2f1f3=interaction(data\$factor_2,data\$factor_1, data\$factor_3)
treatments_f3f1f2=interaction(data\$factor_3,data\$factor_1, data\$factor_2)
data<-data.frame(data,treatments_f1f2, treatments_f1f3, treatments_f2f3, treatments_f1f2f3, treatments_f2f1f3, treatments_f3f1f2)
m4<-aov(response~-1+treatments_f1f2+factor_3*treatments_f1f2,data=data, contrasts=list(treatments_f1f2=contr.sum, factor_3=contr.sum))
m5<-aov(response~-1+treatments_f1f3+factor_2*treatments_f1f3,data=data, contrasts=list(treatments_f1f3=contr.sum, factor_2=contr.sum))
m6<-aov(response~-1+treatments_f2f3+factor_1*treatments_f2f3,data=data, contrasts=list(treatments_f2f3=contr.sum, factor_1=contr.sum))
m7<-aov(response~-1+treatments_f1f2f3,data=data, contrasts=list(treatments_f1f2f3=contr.sum))
m8<-aov(response~-1+treatments_f2f1f3,data=data, contrasts=list(treatments_f2f1f3=contr.sum))
m9<-aov(response~-1+treatments_f3f1f2,data=data, contrasts=list(treatments_f3f1f2=contr.sum))

standard.error<-round(sderrs(m1)[1:nlevels(data\$factor_1)], 4)
factor_1<-levels(data\$factor_1)
rownames(maf1)=NULL
standard.error<-round(sderrs(m2)[1:nlevels(data\$factor_2)],4)
factor_2<-levels(data\$factor_2)
rownames(maf2)=NULL
standard.error<-round(sderrs(m3)[1:nlevels(data\$factor_3)],4)
factor_3<-levels(data\$factor_3)
rownames(maf3)=NULL
standard.error<-round(sderrs(m4)[1:nlevels(data\$treatments_f1f2)],4)
treatments_f1f2<-levels(data\$treatments_f1f2)
rownames(mat1)=NULL
standard.error<-round(sderrs(m5)[1:nlevels(data\$treatments_f1f3)],4)
treatments_f1f3<-levels(data\$treatments_f1f3)
rownames(mat2)=NULL
standard.error<-round(sderrs(m6)[1:nlevels(data\$treatments_f2f3)],4)
treatments_f2f3<-levels(data\$treatments_f2f3)
rownames(mat3)=NULL
standard.error<-round(sderrs(m7)[1:nlevels(data\$treatments_f1f2f3)],4)
treatments_f1f2f3<-levels(data\$treatments_f1f2f3)
rownames(mat4)=NULL
dff=df.residual(m)
test1=fmm(maf1,dff)
test2=fmm(maf2,dff)
test3=fmm(maf3,dff)
groups1=ft(test1, alpha); maf1=maf1[order(maf1[,2], decreasing=TRUE),]
mf1=data.frame(maf1,groups1)
rownames(mf1) = NULL
groups2=ft(test2, alpha); maf2=maf2[order(maf2[,2], decreasing=TRUE),]
mf2=data.frame(maf2,groups2)
rownames(mf2) = NULL
groups3=ft(test3, alpha); maf3=maf3[order(maf3[,2], decreasing=TRUE),]
mf3=data.frame(maf3,groups3)
rownames(mf3) = NULL
c1=rep(1:nlevels(data\$factor_2), each=nlevels(data\$factor_1))
cs1=split(mat1, c1)
test4=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_2));n11=data.frame(n1,levels(data\$factor_2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test4)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2))
cs2=split(mat1, c2)
test5=lapply(cs2, fmm, dff=dff); n2=rep("factor_2 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test5)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups4=lapply(test4, ft, alpha)
mft1=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft1[[x]],groups4[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_2)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups5=lapply(test5, ft, alpha)
mft2=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft2[[x]],groups5[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
c1=rep(1:nlevels(data\$factor_3), each=nlevels(data\$factor_1))
cs1=split(mat2, c1)
test6=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_3));n11=data.frame(n1,levels(data\$factor_3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test6)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_3))
cs2=split(mat2, c2)
test7=lapply(cs2, fmm, dff=dff); n2=rep("factor_3 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test7)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups6=lapply(test6, ft, alpha)
mft3=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft3[[x]],groups6[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_3)))
mft3=lapply(nf1,dd)
names(mft3)=n11
groups7=lapply(test7, ft, alpha)
mft4=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft4[[x]],groups7[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft4=lapply(nf2,ddd)
names(mft4)=n22
c1=rep(1:nlevels(data\$factor_3), each=nlevels(data\$factor_2))
cs1=split(mat3, c1)
test8=lapply(cs1, fmm, dff=dff); n1=rep("factor_2 in", nlevels(data\$factor_3));n11=data.frame(n1,levels(data\$factor_3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test8)=n11
c2=rep(1:nlevels(data\$factor_2), nlevels(data\$factor_3))
cs2=split(mat3, c2)
test9=lapply(cs2, fmm, dff=dff); n2=rep("factor_3 in", nlevels(data\$factor_2));n22=data.frame(n2,levels(data\$factor_2))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test9)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups8=lapply(test8, ft, alpha)
mft5=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft5[[x]],groups8[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_3)))
mft5=lapply(nf1,dd)
names(mft5)=n11
groups9=lapply(test9, ft, alpha)
mft6=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft6[[x]],groups9[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_2)))
mft6=lapply(nf2,ddd)
names(mft6)=n22
c1=rep(1:nlevels(data\$treatments_f2f3), each=nlevels(data\$factor_1))
cs1=split(mat4, c1)
test10=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$treatments_f2f3));n11=data.frame(n1,levels(data\$treatments_f2f3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test10)=n11
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups10=lapply(test10, ft, alpha)
mft7=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft7[[x]],groups10[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$treatments_f2f3)))
mft7=lapply(nf1,dd)
names(mft7)=n11

c1=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2));c1=rep(c1,nlevels(data\$factor_3))
s=rep(0:(nlevels(data\$factor_3)-1), each=nlevels(data\$treatments_f1f2)); oi=max(c1); s=oi*s; c1=c1+s
cs1=split(mat4, c1)
test11=lapply(cs1, fmm, dff=dff); n1=rep("factor_2 in", nlevels(data\$treatments_f1f3));n11=data.frame(n1,levels(data\$treatments_f1f3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test11)=n11
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups11=lapply(test11, ft, alpha)
mft8=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft8[[x]],groups11[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$treatments_f1f3)))
mft8=lapply(nf1,dd)
names(mft8)=n11

c1=rep(1:nlevels(data\$treatments_f1f2), nlevels(data\$factor_3))
cs1=split(mat4, c1)
test12=lapply(cs1, fmm, dff=dff); n1=rep("factor_3 in", nlevels(data\$treatments_f1f2));n11=data.frame(n1,levels(data\$treatments_f1f2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test12)=n11
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups12=lapply(test12, ft, alpha)
mft9=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft9[[x]],groups12[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$treatments_f1f2)))
mft9=lapply(nf1,dd)
names(mft9)=n11
l<-list(a3,  mf1, test1, mf2,test2,mf3,test3,   mft1,test4,mft2,test5,  mft3,test6,mft4,test7,  mft5,test8,mft6,test9,  mft7,test10, mft8,test11, mft9,test12, res)
pres(m)
return(l)}

f8<-function(data, cov){
names(data)=c("factor_1","factor_2", "factor_3", "blocks","response")
data<-data.frame(factor_1=factor(data\$factor_1), factor_2=factor(data\$factor_2), factor_3=factor(data\$factor_3), blocks=factor(data\$blocks), response=data\$response)
m<-aov(response~factor_1*factor_2*factor_3+ blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum, blocks=contr.sum))
m1<-aov(response~-1+ factor_1*factor_2*factor_3+ blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum, blocks=contr.sum))
m2<-aov(response~-1+ factor_2*factor_1*factor_3+ blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum, blocks=contr.sum))
m3<-aov(response~-1+ factor_3*factor_1*factor_2+ blocks,data=data, contrasts=list(factor_1=contr.sum, factor_2=contr.sum, factor_3=contr.sum, blocks=contr.sum))
a2 <- drop1(m,.~.,test="F")
a3 <- a2
a3<-fa2(a3,m)
data2<-na.omit(data)
res=fr3(m,data2)
treatments_f1f2=interaction(data\$factor_1,data\$factor_2)
treatments_f1f3=interaction(data\$factor_1,data\$factor_3)
treatments_f2f3=interaction(data\$factor_2,data\$factor_3)
treatments_f1f2f3=interaction(data\$factor_1,data\$factor_2, data\$factor_3)
treatments_f2f1f3=interaction(data\$factor_2,data\$factor_1, data\$factor_3)
treatments_f3f1f2=interaction(data\$factor_3,data\$factor_1, data\$factor_2)
data<-data.frame(data,treatments_f1f2, treatments_f1f3, treatments_f2f3, treatments_f1f2f3, treatments_f2f1f3, treatments_f3f1f2)
m4<-aov(response~-1+treatments_f1f2+factor_3*treatments_f1f2+ blocks,data=data, contrasts=list(treatments_f1f2=contr.sum, factor_3=contr.sum, blocks=contr.sum))
m5<-aov(response~-1+treatments_f1f3+factor_2*treatments_f1f3+ blocks,data=data, contrasts=list(treatments_f1f3=contr.sum, factor_2=contr.sum, blocks=contr.sum))
m6<-aov(response~-1+treatments_f2f3+factor_1*treatments_f2f3+blocks,data=data, contrasts=list(treatments_f2f3=contr.sum, factor_1=contr.sum, blocks=contr.sum))
m7<-aov(response~-1+treatments_f1f2f3+ blocks,data=data, contrasts=list(treatments_f1f2f3=contr.sum, blocks=contr.sum))
m8<-aov(response~-1+treatments_f2f1f3+ blocks,data=data, contrasts=list(treatments_f2f1f3=contr.sum, blocks=contr.sum))
m9<-aov(response~-1+treatments_f3f1f2+ blocks,data=data, contrasts=list(treatments_f3f1f2=contr.sum, blocks=contr.sum))

standard.error<-round(sderrs(m1) [1:nlevels(data\$factor_1)], 4)
factor_1<-levels(data\$factor_1)
rownames(maf1)=NULL
standard.error<-round(sderrs(m2)[1:nlevels(data\$factor_2)],4)
factor_2<-levels(data\$factor_2)
rownames(maf2)=NULL
standard.error<-round(sderrs(m3)[1:nlevels(data\$factor_3)],4)
factor_3<-levels(data\$factor_3)
rownames(maf3)=NULL
standard.error<-round(sderrs(m4)[1:nlevels(data\$treatments_f1f2)],4)
treatments_f1f2<-levels(data\$treatments_f1f2)
rownames(mat1)=NULL
standard.error<-round(sderrs(m5)[1:nlevels(data\$treatments_f1f3)],4)
treatments_f1f3<-levels(data\$treatments_f1f3)
rownames(mat2)=NULL
standard.error<-round(sderrs(m6)[1:nlevels(data\$treatments_f2f3)],4)
treatments_f2f3<-levels(data\$treatments_f2f3)
rownames(mat3)=NULL
standard.error<-round(sderrs(m7)[1:nlevels(data\$treatments_f1f2f3)],4)
treatments_f1f2f3<-levels(data\$treatments_f1f2f3)
rownames(mat4)=NULL
dff=df.residual(m)
test1=fmm(maf1,dff)
test2=fmm(maf2,dff)
test3=fmm(maf3,dff)
groups1=ft(test1, alpha); maf1=maf1[order(maf1[,2], decreasing=TRUE),]
mf1=data.frame(maf1,groups1)
rownames(mf1) = NULL
groups2=ft(test2, alpha); maf2=maf2[order(maf2[,2], decreasing=TRUE),]
mf2=data.frame(maf2,groups2)
rownames(mf2) = NULL
groups3=ft(test3, alpha); maf3=maf3[order(maf3[,2], decreasing=TRUE),]
mf3=data.frame(maf3,groups3)
rownames(mf3) = NULL
c1=rep(1:nlevels(data\$factor_2), each=nlevels(data\$factor_1))
cs1=split(mat1, c1)
test4=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_2));n11=data.frame(n1,levels(data\$factor_2))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test4)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_2))
cs2=split(mat1, c2)
test5=lapply(cs2, fmm, dff=dff); n2=rep("factor_2 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test5)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups4=lapply(test4, ft, alpha)
mft1=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft1[[x]],groups4[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_2)))
mft1=lapply(nf1,dd)
names(mft1)=n11
groups5=lapply(test5, ft, alpha)
mft2=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft2[[x]],groups5[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft2=lapply(nf2,ddd)
names(mft2)=n22
c1=rep(1:nlevels(data\$factor_3), each=nlevels(data\$factor_1))
cs1=split(mat2, c1)
test6=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$factor_3));n11=data.frame(n1,levels(data\$factor_3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test6)=n11
c2=rep(1:nlevels(data\$factor_1), nlevels(data\$factor_3))
cs2=split(mat2, c2)
test7=lapply(cs2, fmm, dff=dff); n2=rep("factor_3 in", nlevels(data\$factor_1));n22=data.frame(n2,levels(data\$factor_1))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test7)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups6=lapply(test6, ft, alpha)
mft3=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft3[[x]],groups6[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_3)))
mft3=lapply(nf1,dd)
names(mft3)=n11
groups7=lapply(test7, ft, alpha)
mft4=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft4[[x]],groups7[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_1)))
mft4=lapply(nf2,ddd)
names(mft4)=n22
c1=rep(1:nlevels(data\$factor_3), each=nlevels(data\$factor_2))
cs1=split(mat3, c1)
test8=lapply(cs1, fmm, dff=dff); n1=rep("factor_2 in", nlevels(data\$factor_3));n11=data.frame(n1,levels(data\$factor_3))
n11 <- apply(t(n11), 2, paste, collapse = "  ")
names(test8)=n11
c2=rep(1:nlevels(data\$factor_2), nlevels(data\$factor_3))
cs2=split(mat3, c2)
test9=lapply(cs2, fmm, dff=dff); n2=rep("factor_3 in", nlevels(data\$factor_2));n22=data.frame(n2,levels(data\$factor_2))
n22 <- apply(t(n22), 2, paste, collapse = "  ")
names(test9)=n22
csf=function(x){a=x[order(x[,2], decreasing=TRUE),];return(a)}
groups8=lapply(test8, ft, alpha)
mft5=lapply(cs1,csf);
dd<-function(x){l=data.frame(mft5[[x]],groups8[[x]]); return(l)}; nf1=as.list(c(1:nlevels(data\$factor_3)))
mft5=lapply(nf1,dd)
names(mft5)=n11
groups9=lapply(test9, ft, alpha)
mft6=lapply(cs2,csf);
ddd<-function(x){l=data.frame(mft6[[x]],groups9[[x]]); return(l)}; nf2=as.list(c(1:nlevels(data\$factor_2)))
mft6=lapply(nf2,ddd)
names(mft6)=n22
c1=rep(1:nlevels(data\$treatments_f2f3), each=nlevels(data\$factor_1))
cs1=split(mat4, c1)
test10=lapply(cs1, fmm, dff=dff); n1=rep("factor_1 in", nlevels(data\$treatments_f2f3));n11=data.frame(n1,levels(data\$treatments_f2f3))
n11 <- apply(t(n11), 2,