R/cmanova.R

In musto101/wilcox_R: What the package does (short line)

cmanova <-
function(J,K,x,grp=c(1:JK),JK=J*K){
#
# Perform the Choi and Marden
# multivariate one-way rank-based ANOVA
# (Choi and Marden, JASA, 1997, 92, 1581-1590.
#
# x can be a matrix with columns corresponding to groups
# or it can have list mode.
#
# Have a J by K design with J independent levels and K dependent
# measures
#
#
x=elimna(x)
if(is.matrix(x))x<-listm(x)
xx<-list()
nvec<-NA
jk<-0
for(j in 1:J){
for(k in 1:K){
jk<-jk+1
xx[[jk]]<-x[[grp[jk]]]
if(k==1)nvec[j]<-length(xx[[jk]])
}}
N<-sum(nvec)
RVALL<-matrix(0,nrow=N,K)
x<-xx
jk<-0
rmean<-matrix(NA,nrow=J,ncol=K)
for(j in 1:J){
RV<-matrix(0,nrow=nvec[j],ncol=K)
jk<-jk+1
temp1<-matrix(x[[jk]],ncol=1)
for(k in 2:K){
jk<-jk+1
temp1<-cbind(temp1,x[[jk]])
}
X<-temp1
if(j==1)XALL<-X
if(j>1)XALL<-rbind(XALL,X)
n<-nvec[j]
for(i in 1:n){
for (ii in 1:n){
temp3<-sqrt(sum((X[i,]-X[ii,])^2))
if(temp3 != 0)RV[i,]<-RV[i,]+(X[i,]-X[ii,])/temp3
}
RV[i,]<-RV[i,]/nvec[j]
if(j==1 && i==1)sighat<-RV[i,]%*%t(RV[i,])
if(j>1 || i>1)sighat<-sighat+RV[i,]%*%t(RV[i,])
}
}
# Assign ranks to pooled data and compute R bar for each group
for(i in 1:N){
for (ii in 1:N){
temp3<-sqrt(sum((XALL[i,]-XALL[ii,])^2))
if(temp3 != 0)RVALL[i,]<-RVALL[i,]+(XALL[i,]-XALL[ii,])/temp3
}
RVALL[i,]<-RVALL[i,]/N
}
bot<-1-nvec[1]
top<-0
for(j in 1:J){
bot<-bot+nvec[j]
top<-top+nvec[j]
flag<-c(bot:top)
rmean[j,]<-apply(RVALL[flag,],2,mean)
}
sighat<-sighat/(N-J)
shatinv<-solve(sighat)
KW<-0
for(j in 1:J){
KW<-KW+nvec[j]*t(rmean[j,])%*%shatinv%*%rmean[j,]
}
df<-K*(J-1)
sig.level<-1-pchisq(KW,df)
list(test.stat=KW[1,1],df=df,p.value=sig.level)
}