R/sppbb.R

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

sppbb <-
function(J,K,x,est=tmean,JK=J*K,grp=c(1:JK),nboot=500,SEED=TRUE,pr=TRUE,...){
#
# A percentile bootstrap for main effects
# among dependent groups in a split-plot design
# The analysis is done based on all pairs
# of difference scores. The null hypothesis is that
# all such differences have a typical value of zero.
#
#  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(pr)print('As of Oct, 2014, the argument est defaults to tmean')
       if(is.matrix(x)) {
                y <- list()
                for(j in 1:ncol(x))
                        y[[j]] <- x[, j]
                x <- y
}

JK<-J*K
data<-list()
for(j in 1:length(x)){
data[[j]]<-x[[grp[j]]] # Now have the groups in proper order.
}
x<-data
jp<-1-K
kv<-0
kv2<-0
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]
}}
xx<-x
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
# Next determine the n_j values
nvec<-NA
jp<-1-K
for(j in 1:J){
jp<-jp+K
nvec[j]<-length(x[[jp]])
}
#
# Now stack the data in an N by K matrix
#
x<-matrix(NA,nrow=nvec[1],ncol=K)
#
for(k in 1:K)x[,k]<-xx[[k]]
kc<-K
for(j in 2:J){
temp<-matrix(NA,nrow=nvec[j],ncol=K)
for(k in 1:K){
kc<-kc+1
temp[,k]<-xx[[kc]]
}
x<-rbind(x,temp)
}
# Now call function rmdzero to do the analysis
temp<-rmdzero(x,est=est,nboot=nboot,...)
list(p.value=temp$p.value,center=temp$center)
}