How to calculate 95% confidence intervals using bayesboot()

Question

I need help to calculate bootstrap-based credible intervals of the quantity qtt.ci from my regression coef.def.

So far my attempts have resulted in:

Error in quantile.default(s, c(0.025, 0.25, 0.5, 0.75, 0.975)) : missing values and NaN's not allowed if 'na.rm' is FALSE

preceded by:

Warning message: In bayesboot(dat, boot_fn) : The sample from bayesboot contains either NAs, NaNs or NULLs. Make sure that your statistic function only return actual values.

Here are my sample data:

dat <- data.frame(
  A = c(1, 1, 0, 0), B = c(1, 0, 1, 0),
  Pass = c(278, 100, 153, 79), Fail = c(743, 581, 1232, 1731)

Below is my regression. The quantity I want to get the bootstrap-based 95% credible intervals is qtt.ci:

boot_fn <- function(dat) {
           coef.def = unname(coef(glm(cbind(Pass, Fail) ~ A * B, binomial, 
           dat)))
                          }
qtt.ci <- exp(sum(coef.def[2:4])) - exp(coef.def[2]) - exp(coef.def[3]) + 1

Here is my attempt:

bb_ci <- bayesboot(dat, boot_fn)
summary(bb_ci)

Not certain how to get the bootstrap-based confidence intervals for qtt.ci.

Thank you in advance.

EDIT:

Following the answer by @RuiBarradas, I tried doing bootstrap to get the 95% CI for the quantity qtt.ci (which is the quantity for which I want to get the bootstrapped CI), but no success:

library(bayesboot)

boot_fn <- function(dat) {
      coef.def <- unname(coef(glm(cbind(Pass, Fail) ~ A * B, binomial, dat)))
      qtt<- (exp(sum(coef.def[2:4])) - exp(coef.def[2]) - exp(coef.def[3]) + 1) 
      if(all(!is.na(qtt))) qtt else NULL
    }

Runs <- 1e2
qtt.ci <- bayesboot(dat, boot_fn, R = Runs, R2 = Runs)
summary(qtt.ci)

Quantiles:
 statistic    q2.5%     q25%   median     q75%   q97.5%
        V1 2.705878 2.705878 2.705878 2.705878 2.705878

Therefore, this does not give the CI for qtt.ci. The output is simply the point estimate for qtt:

qtt<-(exp(sum(coef.def[2:4])) - exp(coef.def[2]) - exp(coef.def[3]) + 1) 
qtt
[1] 2.705878

Any help would be much appreciated.

Where is the call to `bayesboot()`? – Rui Barradas Dec 16 '18 at 12:51 — Rui Barradas, Dec 16 '18 at 12:51
Just edited the question. Thanks @RuiBarradas. – Krantz Dec 16 '18 at 12:53 — Krantz, Dec 16 '18 at 12:53
Thanks, @MarcoSandri. Just edited. Any solution? – Krantz Dec 16 '18 at 13:47 — Krantz, Dec 16 '18 at 13:47

Rui Barradas · Answer 1 · 2018-12-16T16:23:02.613

2

The following solves the warning issue. I have tested it with much less runs, instead of 4000 just 100.

library(bayesboot)

boot_fn <- function(dat) {
  fit <- glm(cbind(Pass, Fail) ~ A * B, binomial, dat)
  coef.def <- unname(coef(fit))
  if(all(!is.na(coef.def))) coef.def else NULL
}

Runs <- 1e2
bb_ci <- bayesboot(dat, boot_fn, R = Runs, R2 = Runs)
summary(bb_ci)

Edit.

According to the formula in the question and the dialog in comments with the OP, to get the bootstrap-based CI run:

qtt <- exp(sum(bb_ci[2:4])) - exp(bb_ci[2]) - exp(bb_ci[3]) + 1

edited Dec 16 '18 at 16:23

answered Dec 16 '18 at 14:00

Rui Barradas

70,273
8
34
66

Great! Thanks @RuiBarradas. Actually my objective is to get the `bootstrap-based confidence interval` for the quantity `qtt.ci`. I have edited my question to emphasize this and added my attempt to get the `CI` for `qtt.ci` following your answer. Any idea? – Krantz Dec 16 '18 at 14:48
I believe this has the values you want: `as.data.frame(summary(bb_ci))`. – Rui Barradas Dec 16 '18 at 14:59
`V4 = A:B`. That is not `qtt<-(exp(sum(coef.def[2:4])) - exp(coef.def[2]) - exp(coef.def[3]) + 1)`. The quantity for which I want to get the bootstrap-based CI is `qtt`, not `V4`. Any solution? – Krantz Dec 16 '18 at 15:02
@Krantz OK, I think I've got it. In the formula for `qtt` you use `coef.def` but you don't have that value, you assign the return value of the function to `bb_ci`. Just substitute `bb_ci` for `coef-def`. – Rui Barradas Dec 16 '18 at 15:09
Not certain if I get your idea. Could you kindly edit you response? – Krantz Dec 16 '18 at 15:11
@Krantz Is `qtt <- (exp(sum(bb_ci[2:4])) - exp(bb_ci[2]) - exp(bb_ci[3]) + 1)` what you want? – Rui Barradas Dec 16 '18 at 15:52
The quantity for which I want to get the bootstrap-based CI is `qtt <- (exp(sum(bb_ci[2:4])) - exp(bb_ci[2]) - exp(bb_ci[3]) + 1)`, yes. – Krantz Dec 16 '18 at 15:56
Thanks @RuiBarradas. But.ruunning `qtt <- exp(sum(bb_ci[2:4])) - exp(bb_ci[2]) - exp(bb_ci[3]) + 1` gives 3.806e+78, which is not correct. The true CI for `qtt` is around 1.257277 to 4.356140. Any thoughts? – Krantz Dec 16 '18 at 23:00

How to calculate 95% confidence intervals using bayesboot()

1 Answers1