R/rmmcppbtm.R

rmmcppbtm <-
function(x,alpha=.05,con=0,tr=.2,grp=NA,nboot=NA){
#
#   Using the percentile bootstrap method,
#   compute a .95 confidence interval for all linear contasts
#   specified by con, a J by C matrix, where  C is the number of
#   contrasts to be tested, and the columns of con are the
#   contrast coefficients.
#
#   The trimmed means of dependent groups are being compared.
#   By default, 20% trimming is used.
#
#   nboot is the bootstrap sample size. If not specified, a value will
#   be chosen depending on the number of contrasts there are.
#
#   x can be an n by J matrix or it can have list mode
#
#   For alpha=.05, some critical values have been
#   determined via simulations and are used by this function;
#   otherwise an approximation is used.
#
if(!is.list(x) && !is.matrix(x))stop("Data must be stored in a matrix or in list mode.")
if(is.list(x)){
if(is.matrix(con)){
if(length(x)!=nrow(con))stop("The number of rows in con is not equal to the number of groups.")
}}
if(is.list(x)){
# put the data in an n by J matrix
mat<-matrix(0,length(x[[1]]),length(x))
for (j in 1:length(x))mat[,j]<-x[[j]]
}
if(is.matrix(x) && is.matrix(con)){
if(ncol(x)!=nrow(con))stop("The number of rows in con is not equal to the number of groups.")
mat<-x
}
if(is.matrix(x))mat<-x
if(!is.na(sum(grp)))mat<-mat[,grp]
mat<-elimna(mat) # Remove rows with missing values.
J<-ncol(mat)
Jm<-J-1
if(sum(con^2)==0){
d<-(J^2-J)/2
con<-matrix(0,J,d)
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
}}}
d<-ncol(con)
if(is.na(crit) && tr != .2){
print("A critical value must be specified when")
stop("the amount of trimming differs from .2")
}
if(is.na(nboot)){
if(d<=3)nboot<-1000
if(d==6)nboot<-2000
if(d==10)nboot<-4000
if(d==15)nboot<-8000
if(d==21)nboot<-8000
if(d==28)nboot<-10000
}
n<-nrow(mat)
crit<-NA
if(alpha==.05){
if(d==1)crit<-alpha/2
if(d==3){
crit<-.004
if(n>=15)crit<-.006
if(n>=30)crit<-.007
if(n>=40)crit<-.008
if(n>=100)crit<-.009
}
if(d==6){
crit<-.001
if(n>=15)crit<-.002
if(n>=20)crit<-.0025
if(n>=30)crit<-.0035
if(n>=40)crit<-.004
if(n>=60)crit<-.0045
}
if(d==10){
crit<-.00025
if(n>=15)crit<-.00125
if(n>=20)crit<-.0025
}
if(d==15){
crit<-.0005
if(n>=20)crit<-.0010
if(n>=30)crit<-.0011
if(n>=40)crit<-.0016
if(n>=100)crit<-.0019
}
if(d==21){
crit<-.00025
if(n>=20)crit<-.00037
if(n>=30)crit<-.00075
if(n>=40)crit<-.00087
if(n>=60)crit<-.00115
if(n>=100)crit<-.00125
}
if(d==28){
crit<-.0004
if(n>=30)crit<-.0006
if(n>=60)crit<-.0008
if(n>=100)crit<-.001
}
}
if(is.na(crit)){
crit<-alpha/(2*d)
if(n<20)crit<-crit/2
if(n<=10)crit<-crit/2
}
icl<-ceiling(crit*nboot)+1
icu<-ceiling((1-crit)*nboot)
connum<-ncol(con)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
# data is an nboot by n matrix
xbars<-matrix(0,nboot,ncol(mat))
psihat<-matrix(0,connum,nboot)
print("Taking bootstrap samples. Please wait.")
bvec<-bootdep(mat,tr,nboot)
#
# Now have an nboot by J matrix of bootstrap values.
#
test<-1
for (ic in 1:connum){
psihat[ic,]<-apply(bvec,1,bptdpsi,con[,ic])
test[ic]<-sum((psihat[ic,]>0))/nboot
test[ic]<-min(test[ic],1-test[ic])
}
print("Reminder: Test statistic must be less than critical value in order to reject.")
output<-matrix(0,connum,5)
dimnames(output)<-list(NULL,c("con.num","psihat","test","ci.lower","ci.upper"))
tmeans<-apply(mat,2,mean,trim=tr)
psi<-1
for (ic in 1:ncol(con)){
output[ic,2]<-sum(con[,ic]*tmeans)
output[ic,1]<-ic
output[ic,3]<-test[ic]
temp<-sort(psihat[ic,])
output[ic,4]<-temp[icl]
output[ic,5]<-temp[icu]
}
list(output=output,crit=crit,con=con)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.