R/linconm.R

linconm <-
function(x,con=0,est=onestep,alpha=.05,nboot=500,pr=TRUE,...){
#
#   Compute a 1-alpha confidence interval for a set of d linear contrasts
#   involving M-estimators using a bootstrap method. (See Chapter 6.)
#   Independent groups are assumed.
#
#   The data are assumed to be stored in x in list mode.  Thus,
#   x[[1]] contains the data for the first group, x[[2]] the data
#   for the second group, etc. Length(x)=the number of groups = J, say.
#
#   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 measures of location is
#   equal to the sum of the second two, and (2) the difference between
#   the first two is equal to the difference between the measure of location for
#   groups 5 and 6.
#
#   The default number of bootstrap samples is nboot=399
#
#   This function uses the function trimpartt written for this
#   book.
#
#
#
#
if(pr){
print("Note: confidence intervals are adjusted to control FWE")
print("But p-values are not adjusted to control FWE")
}
if(is.matrix(x))x<-listm(x)
con<-as.matrix(con)
if(!is.list(x))stop("Data must be stored in list mode.")
J<-length(x)
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
con[k,id]<-0-1
}}}
if(nrow(con)!=length(x))stop("The number of groups does not match the number of contrast coefficients.")
m1<-matrix(0,J,nboot)
m2<-1 # Initialize m2
mval<-1
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
for(j in 1:J){
mval[j]<-est(x[[j]],...)
xcen<-x[[j]]-est(x[[j]],...)
data<-matrix(sample(xcen,size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
m1[j,]<-apply(data,1,est,...) # A J by nboot matrix.
m2[j]<-var(m1[j,])
}
boot<-matrix(0,ncol(con),nboot)
bot<-1
for (d in 1:ncol(con)){
top<-apply(m1,2,trimpartt,con[,d])
#            A vector of length nboot containing psi hat values
consq<-con[,d]^2
bot[d]<-trimpartt(m2,consq)
boot[d,]<-abs(top)/sqrt(bot[d])
}
testb<-apply(boot,2,max)
ic<-floor((1-alpha)*nboot)
testb<-sort(testb)
psihat<-matrix(0,ncol(con),6)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper","se","p.value"))
for (d in 1:ncol(con)){
psihat[d,1]<-d
psihat[d,2]<-trimpartt(mval,con[,d])
psihat[d,3]<-psihat[d,2]-testb[ic]*sqrt(bot[d])
psihat[d,4]<-psihat[d,2]+testb[ic]*sqrt(bot[d])
psihat[d,5]<-sqrt(bot[d])
pval<-mean((boot[d,]<abs(psihat[d,2])/psihat[d,5]))
psihat[d,6]<-1-pval
}
list(psihat=psihat,crit=testb[ic],con=con)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.