R/yuenbt.R

yuenbt <-
function(x,y,tr=.2,alpha=.05,nboot=599,side=TRUE,nullval=0,pr=TRUE,
plotit=FALSE,op=1,SEED=TRUE){
#
#  Compute a 1-alpha confidence interval for the difference between
#  the trimmed means corresponding to two independent groups.
#  The bootstrap-t method is used.
#
#  The default amount of trimming is tr=.2
#  side=T indicates two-sided method using absolute value of the
#  test statistics within the bootstrap; otherwise the equal-tailed method
#  is used.
#
#  This function uses trimse.
#
side<-as.logical(side)
p.value<-NA
yuenbt<-vector(mode="numeric",length=2)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
x<-x[!is.na(x)]  # Remove missing values in x
y<-y[!is.na(y)]  # Remove missing values in y
xcen<-x-mean(x,tr)
ycen<-y-mean(y,tr)
if(!side){
if(pr)print("NOTE: p-value computed only when side=T")
}
test<-(mean(x,tr)-mean(y,tr))/sqrt(trimse(x,tr=tr)^2+trimse(y,tr=tr)^2)
datax<-matrix(sample(xcen,size=length(x)*nboot,replace=TRUE),nrow=nboot)
datay<-matrix(sample(ycen,size=length(y)*nboot,replace=TRUE),nrow=nboot)
top<-apply(datax,1,mean,tr)-apply(datay,1,mean,tr)
botx<-apply(datax,1,trimse,tr)
boty<-apply(datay,1,trimse,tr)
tval<-top/sqrt(botx^2+boty^2)
if(plotit){
if(op == 1)
akerd(tval)
if(op == 2)
rdplot(tval)
}
if(side)tval<-abs(tval)
tval<-sort(tval)
icrit<-floor((1-alpha)*nboot+.5)
ibot<-floor(alpha*nboot/2+.5)
itop<-floor((1-alpha/2)*nboot+.5)
se<-sqrt((trimse(x,tr))^2+(trimse(y,tr))^2)
yuenbt[1]<-mean(x,tr)-mean(y,tr)-tval[itop]*se
yuenbt[2]<-mean(x,tr)-mean(y,tr)-tval[ibot]*se
if(side){
yuenbt[1]<-mean(x,tr)-mean(y,tr)-tval[icrit]*se
yuenbt[2]<-mean(x,tr)-mean(y,tr)+tval[icrit]*se
p.value<-(sum(abs(test)<=abs(tval)))/nboot
}
list(ci=yuenbt,test.stat=test,p.value=p.value,est.1=mean(x,tr),est.2=mean(y,tr),est.dif=mean(x,tr)-mean(y,tr),
n1=length(x),n2=length(y))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.