R/ydbt.R

ydbt <-
function(x,y,tr=.2,alpha=.05,nboot=599,side=TRUE,plotit=FALSE,op=1,SEED=TRUE){
#
#   Using the bootstrap-t method,
#   compute a .95 confidence interval for the difference between
#   the marginal trimmed means of paired data.
#   By default, 20% trimming is used with B=599 bootstrap samples.
#
#   side=F returns equal-tailed ci
#   side=T returns symmetric ci.
#
side<-as.logical(side)
if(length(x)!=length(y))stop("Must have equal sample sizes.")
m<-cbind(x,y)
m<-elimna(m)
x<-m[,1]
y<-m[,2]
if(sum(c(!is.na(x),!is.na(y)))!=(length(x)+length(y)))stop("Missing values are not allowed.")
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
xcen<-x-mean(x,tr)
ycen<-y-mean(y,tr)
bvec<-apply(data,1,tsub,xcen,ycen,tr)
# bvec is a 1 by nboot matrix containing the bootstrap test statistics.
dotest=yuend(x,y,tr=tr)
estse<-dotest$se
p.value=NULL
dif<-mean(x,tr)-mean(y,tr)
if(!side){
ilow<-round((alpha/2)*nboot)
ihi<-nboot-ilow
bsort<-sort(bvec)
ci<-0
ci[1]<-dif-bsort[ihi]*estse
ci[2]<-dif-bsort[ilow+1]*estse
}
if(side){
bsort<-sort(abs(bvec))
ic<-round((1-alpha)*nboot)
ci<-0
ci[1]<-dif-bsort[ic]*estse
ci[2]<-dif+bsort[ic]*estse
p.value<-(sum(abs(dotest$teststat)<=abs(bvec)))/nboot
}
if(plotit){
if(op==1)akerd(bsort)
if(op==2)rdplot(bsort)
if(op==3)boxplot(bsort)
}
list(ci=ci,dif=dif,p.value=p.value)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.