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.

Rui Barradas Rui Barradas · Accepted Answer · 2018-12-16T14:00:01

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

How to calculate 95% confidence intervals using bayesboot()

1 Answers