Hi,
Many thanks for your responses.
Le samedi 6 novembre 2021, 08:39:22 UTC+1, Rui Barradas <ruipbarra...@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@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.
______________________________________________
R-help@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.
