wilcox_R: What the package does (short line)

t3way <-
function(J,K,L,x,tr=.2,grp=c(1:p),alpha=.05,p=J*K*L,MAT=FALSE,
lev.col=c(1:3),var.col=4,pr=TRUE,IV1=NULL,IV2=NULL,IV3=NULL){
#  Perform a J by K by L (three-way) anova on trimmed means where
#  all JKL groups are independent.
#
#  The R variable data is assumed to contain the raw
#  data stored in list mode. data[[1]] contains the data
#  for the first level of all three factors: level 1,1,1.
#  data][2]] is assumed to contain the data for level 1 of the
#  first two factors and level 2 of the third factor: level 1,1,2
#  data[[L]] is the data for level 1,1,L
#  data[[L+1]] is the data for level 1,2,1. data[[2L]] is level 1,2,L.
#  data[[KL+1]] is level 2,1,1, etc.
#
#  The default amount of trimming is tr=.2
#
#  It is assumed that data has length JKL, the total number of
#  groups being tested.
#
#  MAT=T, assumes data are stored in matrix with 3 columns indicating
#  levels of the three factors.
#  That is, this function calls selby2 for you.
#
if(is.data.frame(x))x=as.matrix(x)
if(!is.null(IV1[1])){
if(is.null(IV2[1]))stop("IV2 is NULL")
if(is.null(IV3[1]))stop("IV3 is NULL")
if(pr)print("Assuming x is a vector containing all of the data; the dependent variable")
xi=elimna(cbind(x,IV1,IV2,IV3))
x=fac2list(xi[,1],xi[,2:4])
J=length(unique(IV1))
K=length(unique(IV2))
L=length(unique(IV3))
p=J*K*L
}
data=x
if(MAT){
if(!is.matrix(data))stop("With MAT=T, data must be a matrix")
if(length(lev.col)!=3)stop("Argument lev.col should have 3 values")
temp=selby2(data,lev.col,var.col)
lev1=length(unique(temp$grpn[,1]))
lev2=length(unique(temp$grpn[,2]))
lev3=length(unique(temp$grpn[,3]))
gv=apply(temp$grpn,2,rank)
gvad=100*gv[,1]+10*gv[,2]+gv[,3]
grp=rank(gvad)
if(pr){
print(paste("Factor 1 has", lev1, "levels"))
print(paste("Factor 2 has", lev2, "levels"))
print(paste("Factor 3 has", lev3, "levels"))
}
if(J!=lev1)warning("J is being reset to the number of levels found")
if(K!=lev2)warning("K is being reset to the number of levels found")
if(L!=lev3)warning("K is being reset to the number of levels found")
J=lev1
K=lev2
L=lev3
data=temp$x
}
if(is.matrix(data))data=listm(data)
if(!is.list(data))stop("Data are not stored in list mode")
if(p!=length(data)){
print("The total number of groups, based on the specified levels, is")
print(p)
print("The number of groups in data is")
print(length(data))
print("Warning: These two values are not equal")
}
tmeans<-0
h<-0
v<-0
for (i in 1:p){
tmeans[i]<-mean(data[[grp[i]]],tr)
h[i]<-length(data[[grp[i]]])-2*floor(tr*length(data[[grp[i]]]))
#    h is the effective sample size
v[i]<-(length(data[[grp[i]]])-1)*winvar(data[[grp[i]]],tr)/(h[i]*(h[i]-1))
#    v contains the squared standard errors
}
v<-diag(v,p,p)   # Put squared standard errors in a diag matrix.
ij<-matrix(c(rep(1,J)),1,J)
ik<-matrix(c(rep(1,K)),1,K)
il<-matrix(c(rep(1,L)),1,L)
jm1<-J-1
cj<-diag(1,jm1,J)
for (i in 1:jm1)cj[i,i+1]<-0-1
km1<-K-1
ck<-diag(1,km1,K)
for (i in 1:km1)ck[i,i+1]<-0-1
lm1<-L-1
cl<-diag(1,lm1,L)
for (i in 1:lm1)cl[i,i+1]<-0-1
alval<-c(1:999)/1000
#  Do test for factor A
cmat<-kron(cj,kron(ik,il))  # Contrast matrix for factor A
Qa<-johan(cmat,tmeans,v,h,alpha)
A.p.value=t3pval(cmat,tmeans,v,h)
# Do test for factor B
cmat<-kron(ij,kron(ck,il))  # Contrast matrix for factor B
Qb<-johan(cmat,tmeans,v,h,alpha)
B.p.value=t3pval(cmat,tmeans,v,h)
# Do test for factor C
cmat<-kron(ij,kron(ik,cl))  # Contrast matrix for factor C
#Qc<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qc<-johan(cmat,tmeans,v,h,alval[i])
if(Qc$teststat>Qc$crit)break
}
C.p.value=irem/1000
# Do test for factor A by B interaction
cmat<-kron(cj,kron(ck,il))  # Contrast matrix for factor A by B
for(i in 1:999){
irem<-i
Qab<-johan(cmat,tmeans,v,h,alval[i])
if(Qab$teststat>Qab$crit)break
}
AB.p.value=irem/1000
# Do test for factor A by C interaction
cmat<-kron(cj,kron(ik,cl))  # Contrast matrix for factor A by C
for(i in 1:999){
irem<-i
Qac<-johan(cmat,tmeans,v,h,alval[i])
if(Qac$teststat>Qac$crit)break
}
AC.p.value=irem/1000
#Qac<-johan(cmat,tmeans,v,h,alpha)
# Do test for factor B by C interaction
cmat<-kron(ij,kron(ck,cl))  # Contrast matrix for factor B by C
#Qbc<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qbc<-johan(cmat,tmeans,v,h,alval[i])
if(Qbc$teststat>Qbc$crit)break
}
BC.p.value=irem/1000
# Do test for factor A by B by C interaction
cmat<-kron(cj,kron(ck,cl))  # Contrast matrix for factor A by B by C
#Qabc<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qabc<-johan(cmat,tmeans,v,h,alval[i])
if(Qabc$teststat>Qabc$crit)break
}
ABC.p.value=irem/1000
list(Qa=Qa$teststat,Qa.crit=Qa$crit,A.p.value=A.p.value,
Qb=Qb$teststat,Qb.crit=Qb$crit,
B.p.value=B.p.value,
Qc=Qc$teststat,Qc.crit=Qc$crit,C.p.value=C.p.value,
Qab=Qab$teststat,Qab.crit=Qab$crit,
AB.p.value=AB.p.value,
Qac=Qac$teststat,Qac.crit=Qac$crit,AC.p.value=AC.p.value,
Qbc=Qbc$teststat,Qbc.crit=Qbc$crit,
BC.p.value=BC.p.value,
Qabc=Qabc$teststat,Qabc.crit=Qabc$crit,ABC.p.value=ABC.p.value)
}