Skip to content

bootstrap confidence intervals

5 messages · David Winsemius, Rui Barradas, varin sacha

varin sacha
# 
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,]
HodgesLehmann(A)
}
boot.out<-boot(data=dat,statistic=med,R=100)
HodgesLehmann(boot.out$t)

boot.ci(boot.out,type="all")
############################
David Winsemius
#
On 11/5/21 1:16 PM, varin sacha via R-help wrote:
# shouldn't this be

HodgesLehmann(temp)  # ???

# makes no sense to extract a bootstrap sample and then return a value calculated on the full dataset
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
Rui Barradas
# 
Hello,


?s 01:36 de 06/11/21, David Winsemius escreveu:
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
2 days later
varin sacha
# 
Hi,

I really thank you a lot for your response.




Le samedi 6 novembre 2021, 02:37:46 UTC+1, David Winsemius <dwinsemius at comcast.net> a ?crit :
On 11/5/21 1:16 PM, varin sacha via R-help wrote:
# shouldn't this be

HodgesLehmann(temp)? # ???

# makes no sense to extract a bootstrap sample and then return a value calculated on the full dataset
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:
?? est lwr.ci upr.ci
449.75 444.25 453.25??? # seems to be cheating
?? est lwr.ci upr.ci
?? 449??? 364??? 500??? # Much closer to the quantiles above
varin sacha
# 
Hi, 

Many thanks for your responses.



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





Hello,


?s 01:36 de 06/11/21, David Winsemius escreveu:
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 at 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.