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:
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
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 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.