R/lindepbt.R

lindepbt <-
function(x, con = NULL, tr = 0.2, alpha = 0.05,nboot=599,dif=TRUE,
SEED=TRUE){
#
# MCP on trimmed means with FWE controlled with Rom's method
# Using a bootstrap-t method.
#
#  dif=T, difference scores are used. And for linear contrasts a simple
#  extension is used.
#
#  dif=F, hypotheses are tested based on the marginal trimmed means.
#
if(SEED)set.seed(2)
if(is.data.frame(x))x=as.matrix(x)
if(is.list(x))x=matl(x)
if(is.null(con))con=con2way(1,ncol(x))$conB # all pairwise
x=elimna(x)
n=nrow(x)
flagcon=F
if(!is.matrix(x))x<-matl(x)
if(!is.matrix(x))stop("Data must be stored in a matrix or in list mode.")
con<-as.matrix(con)
J<-ncol(x)
xbar<-vector("numeric",J)
nval<-nrow(x)
h1<-nrow(x)-2*floor(tr*nrow(x))
df<-h1-1
xbar=apply(x,2,mean,tr=tr)
if(sum(con^2!=0))CC<-ncol(con)
ncon<-CC
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)!=ncol(x))warning("The number of groups does not match the number
 of contrast coefficients.")
ncon<-ncol(con)
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","p.crit","se"))
temp1<-NA
for (d in 1:ncol(con)){
psihat[d,1]<-d
#
#  !dif  Use marginal trimmed means
#
if(!dif){
psihat[d,2]<-sum(con[,d]*xbar)
#
#
sejk<-0
for(j in 1:J){
for(k in 1:J){
djk<-(nval-1)*wincor(x[,j],x[,k], tr)$cov/(h1*(h1-1))
sejk<-sejk+con[j,d]*con[k,d]*djk
}}
sejk<-sqrt(sejk)
test[d,1]<-d
test[d,2]<-sum(con[,d]*xbar)/sejk
test[d,5]<-sejk
#
# now use boostrap-t to determine p-value
#
data<-matrix(sample(n,size=n*nboot,replace=TRUE),nrow=nboot)
xcen=x
for(j in 1:ncol(x))xcen[,j]=xcen[,j]-tmean(x[,j],tr=tr)
bvec=apply(data,1,lindep.sub,xcen,con[,d],tr)
bsort<-sort(abs(bvec))
ic<-round((1-alpha)*nboot)
ci<-0
psihat[d,3]<-psihat[d,2]-bsort[ic]*test[d,5]
psihat[d,4]<-psihat[d,2]+bsort[ic]*test[d,5]
p.value<-mean(abs(test[d,2])<=abs(bvec))
temp1[d]=p.value
}
if(dif){
for(j in 1:J){
if(j==1)dval<-con[j,d]*x[,j]
if(j>1)dval<-dval+con[j,d]*x[,j]
}
temp=trimcibt(dval,tr=tr,alpha=alpha,nboot=nboot,pr=FALSE)
temp1[d]<-temp$p.value #trimci(dval,tr=tr,pr=FALSE)$p.value
test[d,1]<-d
test[d,5]<-trimse(dval,tr=tr)
psihat[d,2]<-mean(dval,tr=tr)
psihat[d,3]<-temp$ci[1] #psihat[,2]-qt(1-test[,4]/2,df)*test[,5]
psihat[d,4]<-temp$ci[2] #psihat[,2]+qt(1-test[,4]/2,df)*test[,5]
}}
#
#   d ends here
#
test[,3]<-temp1
temp2<-order(0-temp1)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2,3]>=zvec)
test[temp2,4]<-zvec
if(flagcon)num.sig<-sum(test[,4]<=test[,5])
if(!flagcon)num.sig<-sum(test[,3]<=test[,4])
list(test=test,psihat=psihat,con=con,num.sig=num.sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.