R/bwimcpES.R

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

bwimcpES <-
function(J,K,x,tr=.2,JK=J*K,grp=c(1:JK),alpha=.05,SEED=TRUE){
#
#  Same as bwimcp only two measures of effect size reported as well.
# First is xi, the measure of effect size returned by yuenv2, which is a robust generalization of
#  Pearson's correlation. The second is a P(difference scores at level A 1 is less than difference scores at level A 2)
#
# Multiple comparisons for interactions
# in a split-plot design.
# The analysis is done by taking difference scores
# among all pairs of dependent groups and
# determining which of
# these differences differ across levels of Factor A
# using trimmed means.
#
#  FWE is controlled via Hochberg's method
# To adjusted p-values, use the function p.adjust
#
# For MOM or M-estimators, use spmcpi which uses a bootstrap method
#
#  The R variable x is assumed to contain the raw
#  data stored in list mode or in a matrix.
#  If 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: level 1,2
#  x[[K]] is the data for level 1,K
#  x[[K+1]] is the data for level 2,1, x[[2K]] is level 2,K, etc.
#
#  If the data are in a matrix, column 1 is assumed to
#  correspond to x[[1]], column 2 to x[[2]], etc.
#
#  When in list mode x is assumed to have length JK, the total number of
#  groups being tested, but a subset of the data can be analyzed
#  using grp
#
       if(is.matrix(x)) {
                y <- list()
                for(j in 1:ncol(x))
                        y[[j]] <- x[, j]
                x <- y
}

JK<-J*K
if(JK!=length(x))stop('Something is wrong. Expected ',JK,' groups but x contains ', length(x), 'groups instead.')
MJ<-(J^2-J)/2
MK<-(K^2-K)/2
JMK<-J*MK
MJMK<-MJ*MK
Jm<-J-1
data<-list()
for(j in 1:length(x)){
data[[j]]<-x[[grp[j]]] # Now have the groups in proper order.
}
x<-data
output<-matrix(0,MJMK,9)
dimnames(output)<-list(NULL,c('A','A','B','B','psihat','p.value','p.crit','EF.xi','EF.WMW'))
jp<-1-K
kv<-0
kv2<-0
test<-NA
for(j in 1:J){
jp<-jp+K
xmat<-matrix(NA,ncol=K,nrow=length(x[[jp]]))
for(k in 1:K){
kv<-kv+1
xmat[,k]<-x[[kv]]
}
xmat<-elimna(xmat)
for(k in 1:K){
kv2<-kv2+1
x[[kv2]]<-xmat[,k]
}}
m<-matrix(c(1:JK),J,K,byrow=T)
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j<jj){
for(k in 1:K){
for(kk in 1:K){
if(k<kk){
ic<-ic+1
output[ic,1]<-j
output[ic,2]<-jj
output[ic,3]<-k
output[ic,4]<-kk
x1<-x[[m[j,k]]]-x[[m[j,kk]]]
x2<-x[[m[jj,k]]]-x[[m[jj,kk]]]
temp<-yuenv2(x1,x2,SEED=SEED)
output[ic,5]<-mean(x1,tr)-mean(x2,tr)
test[ic]<-temp$p.value
output[ic,6]<-test[ic]
output[ic,8]<-temp$Effect.Size
output[ic,9]=cid(x1,x2)$summary.dvals[1]
}}}}}}
ncon<-length(test)
dvec<-alpha/c(1:ncon)
temp2<-order(0-test)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
output[temp2,7]<-zvec
output[,7]<-output[,7]
output
}