bootstrap confidence intervals

Hello,


?s 01:36 de 06/11/21, David Winsemius escreveu:

On 11/5/21 1:16 PM, varin sacha via R-help wrote:

Dear R-experts,

Here is a toy example. How can I get the bootstrap confidence 
intervals working ?

Many thanks for your help

############################
library(DescTools)
library(boot)
A=c(488,437,500,449,364)
dat<-data.frame(A)
med<-function(d,i) {
temp<-d[i,]

# shouldn't this be

HodgesLehmann(temp)? # ???

# makes no sense to extract a bootstrap sample and then return a value 
calculated on the full dataset

HodgesLehmann(A)
}
boot.out<-boot(data=dat,statistic=med,R=100)

I would have imagined that one could simply extract the quantiles of the 
HodgesLehmann at the appropriate tail probabilities:


quantile(boot.out$t, c(0.025, 0.975))
 ??? 2.5%??? 97.5%
400.5000 488.0001


It doesn't seem reasonable to have bootstrap CI's that are much tighter 
than the estimates on the original data:

 > HodgesLehmann(boot.out$t, conf.level=0.95)

 ?? est lwr.ci upr.ci
449.75 444.25 453.25??? # seems to be cheating

 > HodgesLehmann(dat$A, conf.level=0.95)

 ?? est lwr.ci upr.ci
 ?? 449??? 364??? 500??? # Much closer to the quantiles above

This cheating comes from wilcox.test, which is called by HodgesLehman to 
do the calculations. Below is a function calling wilcox.test directly, 
and the bootstrapped intervals are always equal, no matter what way they 
are computed.



A <- c(488, 437, 500, 449, 364)
dat <- data.frame(A)

med <- function(d,i) {
   temp <- d[i, ]
   HodgesLehmann(temp)
}
med2 <- function(d, i, conf.level = 0.95){
   temp <- d[i, ]
   wilcox.test(temp,
               conf.int = TRUE,
               conf.level = Coalesce(conf.level, 0.8),
               exact = FALSE)$estimate
}

set.seed(2021)
boot.out <- boot(data = dat, statistic = med, R = 100)
set.seed(2021)
boot.out2 <- boot(data = dat, statistic = med2, R = 100, conf.level = 0.95)

HodgesLehmann(boot.out$t)
#[1] 452.75
HodgesLehmann(boot.out2$t)
#[1] 452.75

HodgesLehmann(boot.out$t, conf.level = 0.95)
#     est   lwr.ci   upr.ci
#452.7500 447.2500 458.7499
HodgesLehmann(boot.out2$t, conf.level = 0.95)
#     est   lwr.ci   upr.ci
#452.7500 447.2500 458.7499

quantile(boot.out$t, c(0.025, 0.975))
# 2.5% 97.5%
#400.5 494.0
quantile(boot.out2$t, c(0.025, 0.975))
# 2.5% 97.5%
#400.5 494.0

boot.ci(boot.out, type = "all")    # CI's are
boot.ci(boot.out2, type = "all")   # the same



But the bootstrap statistic vectors t are different:



identical(boot.out$t, boot.out2$t)
#[1] FALSE
all.equal(boot.out$t, boot.out2$t)
#[1] "Mean relative difference: 8.93281e-08"




I haven't time to check what is going on in wilcox.test, its source is a 
bit involved, with many if/else statements, maybe I'll come back to this 
but no promises made.


Hope this helps,

Rui Barradas