[R] bootstrap confidence intervals

[varin sacha]

Hi, 

Many thanks for your responses.



Le samedi 6 novembre 2021, 08:39:22 UTC+1, Rui Barradas <ruipbarradas using sapo.pt> a écrit : 





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


______________________________________________
R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.