R/yhbt.R

In musto101/wilcox_R: What the package does (short line)

yhbt <-
function(x,y,tr=.2,alpha=.05,nboot=600,SEED=TRUE,PV=FALSE){
#
#  Compute a 1-alpha confidence interval for the difference between
#  the trimmed means corresponding to two independent groups.
#  The bootstrap-t method with Hall's transformation is used.
#
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)
print("Taking bootstrap samples. Please wait.")
datax<-matrix(sample(xcen,size=length(x)*nboot,replace=TRUE),nrow=nboot)
datay<-matrix(sample(ycen,size=length(y)*nboot,replace=TRUE),nrow=nboot)
val<-NA
for(ib in 1:nboot)val[ib]<-yhall(datax[ib,],datay[ib,],tr=tr,alpha)$test.stat
temp<-yhall(x,y,tr=tr)
sigtil<-temp$sig.tilda
nhat<-temp$nu.tilda
val<-sort(val)
dif<-mean(x,tr=tr)-mean(y,tr=tr)
ilow<-round(alpha*nboot/2)
il<-ilow+1
uval<-nboot-ilow
b.low<-3*((1+nhat*val[il]-nhat/6)^{1/3})/nhat-3/nhat
uval<-nboot-ilow
b.low<-3*((1+nhat*val[il]-nhat/6)^{1/3})/nhat-3/nhat
b.hi<-3*((1+nhat*val[uval]-nhat/6)^{1/3})/nhat-3/nhat
ci.LOW<-dif-sigtil*b.hi
ci.UP<-dif-sigtil*b.low
pv=NULL
if(PV){
#  Determine p-value
pv=1
flag=F
if(dif !=0){
alpha=seq(1:100)/1000
for(i in 1:length(alpha)){
ilow<-round(alpha[i]*nboot/2)
il<-ilow+1
uval<-nboot-ilow
b.low<-3*((1+nhat*val[il]-nhat/6)^{1/3})/nhat-3/nhat
b.hi<-3*((1+nhat*val[uval]-nhat/6)^{1/3})/nhat-3/nhat
ci.low<-dif-sigtil*b.hi
ci.up<-dif-sigtil*b.low
if(ci.low>0 || ci.up<0){
pv=alpha[i]
flag=T
}
if(flag)break
}
if(!flag){
alpha=c(1:99)/100
for(i in 1:length(alpha)){
ilow<-round(alpha[i]*nboot/2)
il<-ilow+1
uval<-nboot-ilow
b.low<-3*((1+nhat*val[il]-nhat/6)^{1/3})/nhat-3/nhat
b.hi<-3*((1+nhat*val[uval]-nhat/6)^{1/3})/nhat-3/nhat
ci.low<-dif-sigtil*b.hi
ci.up<-dif-sigtil*b.low
if(ci.low>0 || ci.up<0){
pv=alpha[i]
flag=T
}
if(flag)break
}
}}}
list(est.dif=dif,conf.interval=c(ci.LOW,ci.UP),p.value=pv)
}