wilcox_R: What the package does (short line)

t2way <-
function(J,K,x,tr=.2,grp=c(1:p),p=J*K,MAT=FALSE,
lev.col=c(1:2),var.col=3,pr=TRUE,IV1=NULL,IV2=NULL){
#  Perform a J by K  (two-way) ANOVA on trimmed means where
#  all groups are independent.
#
#  The R variable x is assumed to contain the raw
#  data stored in list mode, or a matrix with columns
#  corresponding to groups. If stored in list mode, x[[1]] contains the data
#  for the first level of both factors: level 1,1,.
#  x[[2]] is assumed to contain the data for level 1 of the
#  first factor and level 2 of the second factor: level 1,2
#
#  The default amount of trimming is tr=.2
#
#  It is assumed that x has length JK, the total number of
#  groups being tested.
#
#  MAT=T, assumes x are stored in matrix with 3 columns
#  with two of the columns indicated by the argument
#  lev.col
#  specifying the columns of x containing the values of the
#  levels of the two factors.
#  The outcome variable is in column
#  var.col
#  which defaults to column 3
#  That is, this function calls selby2 for you.
#
#  IV1 and IV2: if specified, taken to be the independent variable
#      That is, the group id values
#      and x is assumed to be a vector containing all of the data
#  EXAMPLE: t2way(x=data,IV1=iv1,IV2=iv2)
#  would do a two-way ANOVA based on group id's in iv1 and iv2 and
#  dependent variable data
#
if(is.data.frame(x))data=as.matrix(x)
if(tr==.5){
print("For medians, use med2way if there are no ties")
print("With ties, use linear contrasts in conjunction with medpb")
stop("")
}
if(MAT){
if(!is.matrix(x))stop("With MAT=T, data must be a matrix")
if(length(lev.col)!=2)stop("Argument lev.col should have 3 values")
temp=selby2(x,lev.col,var.col)
lev1=length(unique(temp$grpn[,1]))
lev2=length(unique(temp$grpn[,2]))
gv=apply(temp$grpn,2,rank)
gvad=10*gv[,1]+gv[,2]
grp=rank(gvad)
if(pr){
print(paste("Factor 1 has", lev1, "levels"))
print(paste("Factor 2 has", lev2, "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")
J=lev1
K=lev2
x=temp$x
}
if(!is.null(IV1[1])){
if(is.null(IV2[1]))stop("IV2 is NULL")
if(pr)print("Assuming data is a vector containing all of the data; the dependent variable")
xi=elimna(cbind(x,IV1,IV2))
J=length(unique(xi[,2]))
K=length(unique(xi[,3]))
x=fac2list(xi[,1],xi[,2:3])
}
if(is.matrix(x))x=listm(x)
if(!is.list(x))stop("Data are not stored in list mode")
if(p!=length(x)){
print("The total number of groups, based on the specified levels, is")
print(p)
print("The number of groups is")
print(length(x))
print("Warning: These two values are not equal")
}
tmeans<-0
h<-0
v<-0
for (i in 1:p){
x[[grp[i]]]=elimna(x[[grp[i]]])
tmeans[i]<-mean(x[[grp[i]]],tr)
h[i]<-length(x[[grp[i]]])-2*floor(tr*length(x[[grp[i]]]))
#    h is the effective sample size
if(winvar(x[[grp[i]]],tr)==0)print(paste('The Winsorized variance is zero for group',i))
v[i]<-(length(x[[grp[i]]])-1)*winvar(x[[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)
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
#  Do test for factor A
cmat<-kron(cj,ik)  # Contrast matrix for factor A
alval<-c(1:999)/1000
for(i in 1:999){
irem<-i
Qa<-johan(cmat,tmeans,v,h,alval[i])
if(i==1)dfA=Qa$df
if(Qa$teststat>Qa$crit)break
}
A.p.value=irem/1000
# Do test for factor B
cmat<-kron(ij,ck)  # Contrast matrix for factor B
for(i in 1:999){
irem<-i
Qb<-johan(cmat,tmeans,v,h,alval[i])
if(i==1)dfB=Qb$df
if(Qb$teststat>Qb$crit)break
}
B.p.value=irem/1000
# Do test for factor A by B interaction
cmat<-kron(cj,ck)  # Contrast matrix for factor A by B
for(i in 1:999){
irem<-i
Qab<-johan(cmat,tmeans,v,h,alval[i])
if(i==1)dfAB=Qab$df
if(Qab$teststat>Qab$crit)break
}
AB.p.value=irem/1000
tmeans=matrix(tmeans,J,K,byrow=T)
list(Qa=Qa$teststat,A.p.value=A.p.value, df.A=dfA,
Qb=Qb$teststat,B.p.value=B.p.value,df.B=dfB,
Qab=Qab$teststat,AB.p.value=AB.p.value,df.AB=dfAB,means=tmeans)
}