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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.