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 -- David.> HodgesLehmann(boot.out$t) > > boot.ci(boot.out,type="all") > ############################ > > ______________________________________________ > 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.
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
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:> 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 -- David.> HodgesLehmann(boot.out$t) > > boot.ci(boot.out,type="all")> ############################ > > ______________________________________________ > 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.______________________________________________ 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.