wilcox_R: What the package does (short line)

bwrank <-
function(J,K,x,grp=c(1:p),p=J*K){
#
# Between by within rank-based ANOVA
# That is, have a J by K design with J independent levels and K dependent
# measures
#
# x can be a matrix with columns corresponding to groups
# or it can have list mode.
#
#
if(is.data.frame(x))data=as.matrix(x)
if(is.matrix(x))x<-listm(x)
x=x[grp]
xx<-list()
nvec<-NA
alldat<-NA
klow<-1-K
kup<-0
iall=0
for (j in 1:J){
klow<-klow+K
kup<-kup+K
mtemp=elimna(matl(x[klow:kup]))
for(k in 1:K){
iall=iall+1
xx[[iall]]=mtemp[,k]
}}
for(j in 1:p){
alldat<-c(alldat,xx[[j]])
nvec[j]<-length(xx[[j]])
}
#
#  Check sample sizes
#
nmat<-matrix(nvec,J,K,byrow=T)
for(j in 1:J){
if(var(nmat[j,]) !=0){
warning("Number of observations among dependent groups for level",paste(j)," of Factor A are unequal")
print("Matrix of sample sizes:")
print(nmat)
}}
if(sum(is.na(alldat[2:length(alldat)])>0))stop("Missing values not allowed")
rval<-rank(alldat[2:length(alldat)])
rdd<-mean(rval) # R bar ...
xr<-list()
il<-1-nvec[1]
iu<-0
for(j in 1:p){
il<-il+nvec[j]
iu<-iu+nvec[j]
xr[[j]]<-rval[il:iu]
}
v<-matrix(0,p,p)
Ja<-matrix(1,J,J)
Ia<-diag(1,J)
Pa<-Ia-Ja/J
Jb<-matrix(1,K,K)
Ib<-diag(1,K)
Pb<-Ib-Jb/K
cona<-kron(Pa,Ib)
conb<-kron(Ia,Pb)
conab<-kron(Pa,Pb)
for(k in 1:K){
temp<-x[[k]]
bigm<-matrix(temp,ncol=1)
jk<-k
for (j in 2:J){
jk<-jk+K
tempc<-matrix(x[[jk]],ncol=1)
bigm<-rbind(bigm,tempc)
temp<-c(temp,x[[jk]])
}}
N<-length(temp)
rbbd<-NA
for(k in 1:K){
bigm<-xr[[k]]
jk<-k
for (j in 2:J){
jk<-jk+K
bigm<-c(bigm,xr[[jk]])
}}
rbjk<-matrix(NA,nrow=J,ncol=K) #R_.jk
ic<-0
for (j in 1:J){
for(k in 1:K){
ic<-ic+1
rbjk[j,k]<-mean(xr[[ic]]) # R bar_.jk
}}
for(k in 1:K)rbbd[k]<-mean(rbjk[,k])
rbj<-1 # R_.j.
sigv<-0
njsam<-0 # n_j
icc<-1-K
ivec<-c(1:K)-K
for (j in 1:J){
icc<-icc+K
ivec<-ivec+K
temp<-xr[ivec[1]:ivec[K]]
temp<-matl(temp)
tempv<-apply(temp,1,mean)
njsam[j]<-nvec[icc]
rbj[j]<-mean(rbjk[j,])
sigv[j]<-var(tempv) # var of R bar_ij.
}
nv<-sum(njsam)
phat<-(rbjk-.5)/(nv*K)
sv2<-sum(sigv/njsam)
uv<-sum((sigv/njsam)^2)
dv<-sum((sigv/njsam)^2/(njsam-1))
testA<-J*var(rbj)/sv2
klow<-1-K
kup<-0
sv<-matrix(0,nrow=K,ncol=K)
rvk<-NA
for(j in 1:J){
klow<-klow+K
kup<-kup+K
sel<-c(klow:kup)
m<-matl(xr[klow:kup])
m<-elimna(m)
xx<-listm(m)
xx<-listm(m)
vsub<-nv*var(m)/(nv*K*nv*K*njsam[j])
v[sel,sel]<-vsub
sv<-sv+vsub
}
sv<-sv/J^2
testB<-nv/(nv*K*nv*K*sum(diag(Pb%*%sv)))*sum((rbbd-mean(rbbd))^2)
testAB<-0
for (j in 1:J){
for (k in 1:K){
testAB<-testAB+(rbjk[j,k]-rbj[j]-rbbd[k]+rdd)^2
}}
testAB<-nv*testAB/(nv*K*nv*K*sum(diag(conab%*%v)))
nu1B<-(sum(diag(Pb%*%sv)))^2/sum((diag(Pb%*%sv%*%Pb%*%sv)))
nu1A<-(J-1)^2/(1+J*(J-2)*uv/sv2^2)
nu2A<-sv2^2/dv
nu1AB<-(sum(diag(conab%*%v)))^2/sum(diag(conab%*%v%*%conab%*%v))
sig.A<-1-pf(testA,nu1A,nu2A)
sig.B<-1-pf(testB,nu1B,1000000)
sig.AB<-1-pf(testAB,nu1AB,1000000)
list(test.A=testA,sig.A=sig.A,test.B=testB,sig.B=sig.B,test.AB=testAB,
sig.AB=sig.AB,avg.ranks=rbjk,rel.effects=phat)
}