wilcox_R: What the package does (short line)

MULtr.anova <-
function(x,J=NULL,p=NULL,tr=.2,alpha=.05){
#
#  Do Multivariate ANOVA with trimmed means using
#  Johansen's method
#
#  x is assumed to have list mode with length(x)=J=number of groups and
#  x[[j]] is an n_j-by-p  matrix, p is the number of variables.
#
#  x can also be a matrix when J and p are specified. It is assumed the data are stored in
#   a matrix in the same manner expected by bwtrim.
#
#  To get a p-value, use the function MULAOVp
#
if(is.matrix(x) || is.data.frame(x)){
if(is.null(J) && is.null(p))stop("Specify J or P")
x=MAT2list(x,p=p,J=J)
}
x=lapply(x,as.matrix)
x=lapply(x,elimna)
p=ncol(x[[1]])
iden=diag(p)
J=length(x)
tvec=list()
nval=lapply(x,nrow)
Rtil=lapply(x,wincov,tr=tr)
tvec=lapply(x,mmean,tr=tr)
g=list()
gmean=rep(0,p) # grand mean eventually
groupm=list()
Wsum=matrix(0,ncol=p,nrow=p)
W=list()
f=0
Aw=0
for(j in 1:J){
tvec[[j]]=as.matrix(tvec[[j]])
g[[j]]=floor(nval[[j]]*tr)
Rtil[[j]]=Rtil[[j]]*(nval[[j]]-1)/((nval[[j]]-2*g[[j]])*(nval[[j]]-2*g[[j]]-1))
f[j]=nval[[j]]-2*g[[j]]-1
W[[j]]=solve(Rtil[[j]])
groupm[[j]]=apply(x[[j]],2,tmean,tr=tr)
Wsum=Wsum+W[[j]]
gmean=gmean+W[[j]]%*%tvec[[j]]
}
Wsuminv=solve(Wsum)
for(j in 1:J){
temp=iden-Wsuminv%*%W[[j]]
tempsq=temp%*%temp
Aw=Aw+(sum(diag(tempsq))+(sum(diag(temp)))^2)/f[j]
}
Aw=Aw/2
gmean=as.matrix(gmean)
gmean=solve(Wsum)%*%gmean # Final weighted grand mean
df=p*(J-1)
crit<-qchisq(1-alpha,df)
crit<-crit+(crit/(2*df))*(Aw+3*Aw*crit/(df+2))
test=0
for(k in 1:p){
for(m in 1:p){
for(j in 1:J){
test=test+W[[j]][k,m]*(groupm[[j]][m]-gmean[m])*(groupm[[j]][k]-gmean[k])
}}}
list(test.stat=test,crit.value=crit)
}