R/lincdm.R

lincdm <-
function(x,con=0,alpha=.05,q=.5,mop=FALSE,nboot=100,SEED=TRUE){
#
#  A heteroscedastic test of d linear contrasts among
#  dependent groups using medians.
#
#  The data are assumed to be stored in $x$ in list mode.
#  Length(x) is assumed to correspond to the total number of groups, J
#  It is assumed all groups are independent.
#
#  con is a J by d matrix containing the contrast coefficients that are used.
#  If con is not specified, all pairwise comparisons are made.
#
#  q is the quantile used to compare groups.
#  con contains contrast coefficients,
#  con=0 means all pairwise comparisons are used
#  mop=F, use single order statistic
#  mop=T, use usual sample median, even if q is not equal to .5
#  in conjunction with a bootstrap estimate of covariances among
#  the medians using
#  nboot samples.
#
#  Missing values are automatically removed.
#
#
if(mop && SEED)set.seed(2)
if(is.list(x)){
x<-matl(x)
x<-elimna(x)
}
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in a matrix or in list mode.")
con<-as.matrix(con)
J<-length(x)
h<-length(x[[1]])
w<-vector("numeric",J)
xbar<-vector("numeric",J)
for(j in 1:J){
if(!mop)xbar[j]<-qest(x[[j]],q=q)
if(mop)xbar[j]<-median(x[[j]])
}
if(sum(con^2)==0){
temp<-qdmcp(x,alpha=alpha,q=q,pr=FALSE)
test<-temp$test
psihat<-temp$psihat
num.sig<-temp$num.sig
}
if(sum(con^2)>0){
ncon<-ncol(con)
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(nrow(con)!=length(x)){
stop("The number of groups does not match the number of contrast coefficients.")
}
psihat<-matrix(0,ncol(con),4)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper"))
test<-matrix(0,ncol(con),5)
dimnames(test)<-list(NULL,c("con.num","test","p.value","crit.p.value","se"))
df<-length(x[[1]])-1
if(!mop)w<-covmmed(x,q=q)
if(mop)w<-bootcov(x,nboot=nboot,pr=FALSE)
for (d in 1:ncol(con)){
psihat[d,1]<-d
psihat[d,2]<-sum(con[,d]*xbar)
cvec<-as.matrix(con[,d])
sejk<-sqrt(t(cvec)%*%w%*%cvec)
test[d,1]<-d
test[d,2]<-sum(con[,d]*xbar)/sejk
test[d,3]<-2*(1-pt(abs(test[d,2]),df))
test[d,5]<-sejk
}
temp1<-test[,3]
temp2<-order(0-temp1)
zvec<-dvec[1:ncon]
test[temp2,4]<-zvec
psihat[,3]<-psihat[,2]-qt(1-test[,4]/2,df)*test[,5]
psihat[,4]<-psihat[,2]+qt(1-test[,4]/2,df)*test[,5]
num.sig<-sum(test[,3]<=test[,4])
}
list(test=test,psihat=psihat,num.sig=num.sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.