R/linconMpb.R

linconMpb <-
function(x,alpha=.05,nboot=1000,grp=NA,est=tmean,con=0,bhop=FALSE,
SEED=TRUE,PDIS=FALSE,J=NULL,p=NULL,...){
#
#   Multiple comparisons for  J independent groups using trimmed means
#   with multivariate data for each group.
#
#   A percentile bootstrap method with Rom's method is used.
#
#   The data are assumed to be stored in x
#   which  has list mode,
#   x[[1]] contains the data for the first group in the form of a
#   matrix, x[[2]] the data
#   for the second group, etc. Length(x)=the number of groups = J.
#
#   est is the measure of location and defaults to the median
#   ... can be used to set optional arguments associated with est
#
#   The argument grp can be used to analyze a subset of the groups
#   Example: grp=c(1,3,5) would compare groups 1, 3 and 5.
#
#   Missing values are automatically removed.
#
con<-as.matrix(con)
if(is.matrix(x) || is.data.frame(x)){
if(is.null(J) && is.null(p))stop("Specify J or P")
x=MAT2list(x,p=p,J=J)
}
if(!is.list(x))stop("Data must be stored in list mode.")
if(!is.na(sum(grp))){  # Only analyze specified groups.
xx<-list()
for(i in 1:length(grp))xx[[i]]<-x[[grp[i]]]
x<-xx
}
J<-length(x)
nullvec=rep(0,ncol(x[[1]]))
bplus=nboot+1
tempn<-0
mvec<-list
for(j in 1:J){
x[[j]]<-elimna(x[[j]])
}
Jm<-J-1
#
# Determine contrast matrix
#
if(sum(con^2)==0){
ncon<-(J^2-J)/2
con<-matrix(0,J,ncon)
id<-0
for (j in 1:Jm){
jp<-j+1
for (k in jp:J){
id<-id+1
con[j,id]<-1
con[k,id]<-0-1
}}}
ncon<-ncol(con)
if(nrow(con)!=J)stop("Something is wrong with con; the number of rows does not match the number of groups.")
# Determine critical levels
if(!bhop){
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(bhop)dvec<-(ncon-c(1:ncon)+1)*alpha/ncon
bvec<-array(NA,c(J,nboot,ncol(x[[1]])))
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
#print("Taking bootstrap samples. Please wait.")
nvec=lapply(x,nrow)
for(j in 1:J){
data<-matrix(sample(nvec[[j]],size=nvec[[j]]*nboot,replace=TRUE),nrow=nboot)
bvec[j,,]<-apply(data,1,linconMpb.sub,x[[j]],est,...) # Bootstrapped values for jth group
}
test<-NA
for (d in 1:ncon){
tv=matrix(0,nboot,ncol(x[[1]])) #nboot by p matrix reflecting Psi hat
estit=rep(0,ncol(x[[1]]))
for(j in 1:J){
tv=tv+con[j,d]*bvec[j,,]
estit=estit+con[j,d]*apply(x[[j]],2,est,...)
}
if(!PDIS)m1=cov(tv)
tv=rbind(tv,nullvec)
if(!PDIS)dv=mahalanobis(tv,center=estit,m1)
if(PDIS)dv=pdis(tv,center=estit) # projection distances
test[d]=1-sum(dv[bplus]>=dv[1:nboot])/nboot
}
output<-matrix(0,ncon,3)
dimnames(output)<-list(NULL,c("con.num","p.value","p.crit"))
temp2<-order(0-test)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
output[temp2,3]<-zvec
for (ic in 1:ncol(con)){
output[ic,1]<-ic
output[ic,2]<-test[ic]
}
num.sig<-sum(output[,2]<=output[,3])
list(output=output,con=con,num.sig=num.sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.