R/linpbg.R

linpbg <-
function(x,con=0,alpha=.05,nboot=NA,est=mest,...){
#
#   Compute a 1-alpha confidence interval
#   for a set of d linear contrasts
#   involving trimmed means using the percentile bootstrap method.
#   Independent groups are assumed.
#
#   The data are assumed to be stored in x in list mode or in a matrix.
#   Thus,
#   x[[1]] contains the data for the first group, x[[2]] the data
#   for the second group, etc.
#   If x has list mode, length(x)=the number of groups = J, say.
#
#   Missing values are automatically removed.
#
#   con is a J by d matrix containing the
#   contrast coefficents of interest.
#   If unspecified, all pairwise comparisons are performed.
#   For example, con[,1]=c(1,1,-1,-1,0,0)
#   and con[,2]=c(,1,-1,0,0,1,-1)
#   will test two contrasts: (1) the sum of the first
#   two trimmed means is
#   equal to the sum of the second two,
#   and (2) the difference between
#   the first two is equal to the difference
#   between the trimmed means of
#   groups 5 and 6.
#
#
con<-as.matrix(con)
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in a matrix or in list mode.")
J<-length(x)
for(j in 1:J){
xx<-x[[j]]
xx[[j]]<-xx[!is.na(xx)] # Remove any missing values.
}
Jm<-J-1
d<-(J^2-J)/2
if(sum(con^2)==0){
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 #If con not specified do all pairwise comparisons
con[k,id]<-0-1
}}}
if(nrow(con)!=length(x)){
stop("The number of groups does not match the number of contrast coefficients.")
}
if(is.na(nboot)){
nboot<-5000
if(ncol(con)<=4)nboot<-2000
}
m1<-matrix(0,nrow=J,ncol=nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
for(j in 1:J){
paste("Working on group ",j)
data<-matrix(sample(x[[j]],size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
m1[j,]<-apply(data,1,est,...)
}
testb<-NA
boot<-matrix(0,ncol(con),nboot)
testvec<-NA
for (d in 1:ncol(con)){
boot[d,]<-apply(m1,2,trimpartt,con[,d])
# A vector of length nboot containing psi hat values
# and corresponding to the dth linear contrast
testb[d]<-sum((boot[d,]>0))/nboot
testvec[d]<-min(testb[d],1-testb[d])
}
#
#  Determine critical value
#
dd<-ncol(con)
if(alpha==.05){
if(dd==1)crit<-alpha/2
if(dd==2)crit<-.014
if(dd==3)crit<-.0085
if(dd==4)crit<-.007
if(dd==5)crit<-.006
if(dd==6)crit<-.0045
if(dd==10)crit<-.0023
if(dd==15)crit<-.0016
}
else{
crit<-alpha/(2*dd)
}
icl<-round(crit*nboot)
icu<-round((1-crit)*nboot)
psihat<-matrix(0,ncol(con),4)
test<-matrix(0,ncol(con),3)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper"))
dimnames(test)<-list(NULL,c("con.num","test","crit.val"))
for (d in 1:ncol(con)){
test[d,1]<-d
psihat[d,1]<-d
testit<-lincon(x,con[,d],tr)
test[d,2]<-testvec[d]
temp<-sort(boot[d,])
psihat[d,3]<-temp[icl]
psihat[d,4]<-temp[icu]
psihat[d,2]<-testit$psihat[1,2]
test[d,3]<-crit
}
list(psihat=psihat,test=test,con=con)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.