R/trimpb2.R

trimpb2 <-
function(x,y,tr=.2,alpha=.05,nboot=2000,WIN=FALSE,win=.1,plotit=FALSE,op=4,
SEED=TRUE){
#
#   Compute a 1-alpha confidence interval for
#   the difference between two 20% trimmed means.
#   Independent groups are assumed.
#
#   The default number of bootstrap samples is nboot=2000
#
#   tr is the amount of trimming
#
#   win is the amount of Winsorizing before bootstrapping
#   when WIN=T.
#
#   Missing values are automatically removed.
#
x<-x[!is.na(x)]
y<-y[!is.na(y)]
if(WIN){
if(win>tr)stop("Cannot Winsorize more than you trim")
if(tr < .2){print("When Winsorizing, the amount of trimming")
print("should be at least .2")
}
if(min(c(length(x),length(y))) < 15){
print ("Warning: Winsorizing with sample sizes less than 15")
print("can result in poor control over the probability of a Type I error")
}
x<-winval(x,win)
y<-winval(y,win)
}
xx<-list()
xx[[1]]<-x
xx[[2]]<-y
est.dif<-tmean(xx[[1]],tr=tr)-tmean(xx[[2]],tr=tr)
crit<-alpha/2
temp<-round(crit*nboot)
icl<-temp+1
icu<-nboot-temp
bvec<-matrix(NA,nrow=2,ncol=nboot)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
for(j in 1:2){
data<-matrix(sample(xx[[j]],size=length(xx[[j]])*nboot,replace=TRUE),nrow=nboot)
bvec[j,]<-apply(data,1,mean,tr) # Bootstrapped trimmed means for jth group
}
top<-bvec[1,]-bvec[2,]
test<-sum(top<0)/nboot+.5*sum(top==0)/nboot
if(test > .5)test<-1-test
top<-sort(top)
ci<-NA
ci[1]<-top[icl]
ci[2]<-top[icu]
if(plotit)g2plot(bvec[1,],bvec[2,],op=op)
list(p.value=2*test,ci=ci,est.dif=est.dif)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.