R/mcp2a.R

mcp2a <-
function(J,K,x,est=mom,con=NULL,alpha=.05,nboot=NA,grp=NA,...){
#
# Do all pairwise comparisons of
# main effects for Factor A and B and all interactions
#
        #  The data are assumed to be stored in x
        #  in list mode or in a matrix.
        #  If grp is unspecified, it is assumed 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 factor: level 1,2
        #  x[[j+1]] is the data for level 2,1, etc.
        #  If the data are in wrong order, grp can be used to rearrange the
        #  groups. For example, for a two by two design, grp<-c(2,4,3,1)
        #  indicates that the second group corresponds to level 1,1;
        #  group 4 corresponds to level 1,2; group 3 is level 2,1;
        #  and group 1 is level 2,2.
        #
        #   Missing values are automatically removed.
        #
        JK <- J * K
        if(is.matrix(x))
                x <- listm(x)
        if(!is.na(grp)) {
                yy <- x
                for(j in 1:length(grp))
                        x[[j]] <- yy[[grp[j]]]
        }
        if(!is.list(x))
                stop("Data must be stored in list mode or a matrix.")
mvec<-NA
  tempn=0
        for(j in 1:JK) {
                xx <- x[[j]]
                x[[j]] <- xx[!is.na(xx)]
                mvec[j]<-est(x[[j]],...)
tempn[j]=length(x[[j]])
        }
nmax=max(tempn)
        #
        # Create the three contrast matrices
        #
        if(JK != length(x))
                warning("The number of groups does not match the number of contrast coefficients.")
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
#  Determine nboot if a value was not specified
if(is.na(nboot)){
nboot<-5000
if(J <= 8)nboot<-4000
if(J <= 3)nboot<-2000
}
bvec<-matrix(NA,nrow=JK,ncol=nboot)
for(j in 1:JK){
data<-matrix(sample(x[[j]],size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
bvec[j,]<-apply(data,1,est,...) # J by nboot matrix, jth row contains
#                          bootstrapped  estimates for jth group
}
outvec<-list()
if(!is.null(con))stop('Use linconm when specifying the linear contrast coefficients')
temp3<-con2way(J,K)
for(jj in 1:3){
con<-temp3[[jj]]
con<-as.matrix(con)
ncon<-ncol(con)
# Determine critical values
if(alpha==.05){
dvec<-c(.025,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.005,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(nmax>80){
if(alpha==.05){
dvec<-c(.05,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.01,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01){
dvec<-alpha/c(1:ncon)
}}
if(alpha != .05 && alpha != .01)dvec<-alpha/c(1:ncon)
test<-NA
bcon<-t(con)%*%bvec #ncon by nboot matrix
tvec<-t(con)%*%mvec
for (d in 1:ncon){
test[d]<-sum(bcon[d,]>0)/nboot
if(test[d]> .5)test[d]<-1-test[d]
}
output<-matrix(0,ncon,6)
dimnames(output)<-list(NULL,c("con.num","psihat","sig.test","sig.crit","ci.lower","ci.upper"))
temp2<-order(0-test)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
if(sum(sigvec)<ncon){
dd<-ncon-sum(sigvec) #number that are sig.
ddd<-sum(sigvec)+1
zvec[ddd:ncon]<-dvec[ddd]
}
output[temp2,4]<-zvec
icl<-round(dvec[ncon]*nboot)+1
icu<-nboot-icl-1
for (ic in 1:ncol(con)){
output[ic,2]<-tvec[ic,]
output[ic,1]<-ic
output[ic,3]<-test[ic]
temp<-sort(bcon[ic,])
output[ic,5]<-temp[icl]
output[ic,6]<-temp[icu]
}
outvec[[jj]]<-output
}
list(FactorA=outvec[[1]],FactorB=outvec[[2]],Interactions=outvec[[3]],
conA=temp3[[1]],conB=temp3[[2]],conAB=temp3[[3]])
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.