Bug? A minimal working example for bcajack2(Blist)

[t-pollington]

Dear @bnaras,

I'm having trouble writing an MWE for bcajack2(). This was for my own benefit but also was planning to do a pull request to update the bcajack2() example as it doesn't really show off its specific functionality i.e. being able to do separate \hat{theta}* bootstrap calculations elsewhere and then get a 95% BCa CI, which is the use case I am most interested in. I'm in the middle of writing paper where I hope to do just this and would love to use your bcaboot package and cite it too!

From reading the function code I believe the error emanates from the internal function qbca2() which provides vl0$lims. On line 98 log(Yi) will produce -Inf when Yi[i] = 0; could this be the cause? Since I noticed when not using the Blist format on line 165 that Y gets values from table() which would not include values not present within a range i.e. table(c(1,2,4)) doesn't include a zero for 3 and therefore later doesn't throw the same error as Blist sees?

data(diabetes, package = "bcaboot")
Xy <- cbind(diabetes$x, diabetes$y)
rfun <- function(Xy) {
  y <- Xy[, 11]
  X <- Xy[, 1:10]
  summary(lm(y~X) )$adj.r.squared
}

# letting bcajack2 perform the bootstrapping
set.seed(1234)
bcaboot::bcajack2(x = Xy, B = 1000, func = rfun, m = 40, verbose = FALSE)
#> $call
#> bcaboot::bcajack2(x = Xy, B = 1000, func = rfun, m = 40, verbose = FALSE)
#> 
#> $lims
#>             bca      jacksd       std   pct
#> 0.025 0.4289798 0.003371377 0.4462564 0.002
#> 0.05  0.4374506 0.006440959 0.4559517 0.006
#> 0.1   0.4533929 0.004038199 0.4671297 0.018
#> 0.16  0.4621266 0.003321927 0.4759630 0.037
#> 0.5   0.4950141 0.001790731 0.5065603 0.219
#> 0.84  0.5255210 0.002355254 0.5371577 0.586
#> 0.9   0.5357038 0.002870023 0.5459909 0.692
#> 0.95  0.5459683 0.002031757 0.5571690 0.804
#> 0.975 0.5540292 0.003463861 0.5668642 0.877
#> 
#> $stats
#>         theta       sdboot          z0             a      sdjack
#> est 0.5065603 0.0307678772 -0.38532047 -0.0109424195 0.031399281
#> jsd 0.0000000 0.0007431751  0.03769735  0.0007008304 0.000282446
#> 
#> $B.mean
#> [1] 1000.0000000    0.5189562
#> 
#> $ustats
#>      ustat        sdu 
#> 0.49416441 0.03171256 
#> 
#> attr(,"class")
#> [1] "bcaboot"

# giving backjack2 the precomputed bootstrap estimates----
set.seed(1234)
inds = matrix(NA, nrow = 1000, ncol = nrow(Xy))
inds = t(apply(inds, MARGIN = 1, function(inds) sample(1:nrow(Xy), replace = TRUE))) # get bootstrap indices to feed into estimator rfun()
tt = vector(mode = "numeric", length = 1000)
tt = apply(inds, MARGIN = 1, function(inds) rfun(Xy[inds,])) # get \hat{theta}*
Y = t(apply(inds, MARGIN = 1, tabulate, nbins = nrow(Xy)))

Blist = list(Y = Y, 
             tt = tt, 
             t0 = rfun(Xy)) # original estimate \hat{theta}
bcaboot::bcajack2(Blist) # "argument "B" is missing, with no default" which contradicts https://github.com/bnaras/bcaboot/blob/master/man/bcajack2.Rd line 114
#> Error in bcaboot::bcajack2(Blist): argument "B" is missing, with no default
bcaboot::bcajack2(B = Blist) # bca $lims are all NA
#> $call
#> bcaboot::bcajack2(B = Blist)
#> 
#> $lims
#>       bca jacksd       std pct
#> 0.025  NA    NaN 0.4430264  NA
#> 0.05   NA    NaN 0.4532409  NA
#> 0.1    NA    NaN 0.4650177  NA
#> 0.16   NA    NaN 0.4743241  NA
#> 0.5    NA    NaN 0.5065603  NA
#> 0.84   NA    NaN 0.5387965  NA
#> 0.9    NA    NaN 0.5481029  NA
#> 0.95   NA    NaN 0.5598797  NA
#> 0.975  NA    NaN 0.5700943  NA
#> 
#> $stats
#>         theta       sdboot          z0   a sdjack
#> est 0.5065603 0.0324158820 -0.20701262  NA     NA
#> jsd 0.0000000 0.0007799418  0.03852254 NaN    NaN
#> 
#> $B.mean
#> [1] 1000.0000000    0.5137792
#> 
#> $ustats
#>     ustat       sdu 
#> 0.4993414        NA 
#> 
#> attr(,"class")
#> [1] "bcaboot"
bcaboot::bcajack2(B = Blist, m = 40) # still NA
#> $call
#> bcaboot::bcajack2(B = Blist, m = 40)
#> 
#> $lims
#>       bca jacksd       std pct
#> 0.025  NA    NaN 0.4430264  NA
#> 0.05   NA    NaN 0.4532409  NA
#> 0.1    NA    NaN 0.4650177  NA
#> 0.16   NA    NaN 0.4743241  NA
#> 0.5    NA    NaN 0.5065603  NA
#> 0.84   NA    NaN 0.5387965  NA
#> 0.9    NA    NaN 0.5481029  NA
#> 0.95   NA    NaN 0.5598797  NA
#> 0.975  NA    NaN 0.5700943  NA
#> 
#> $stats
#>         theta      sdboot          z0   a sdjack
#> est 0.5065603 0.032415882 -0.20701262  NA     NA
#> jsd 0.0000000 0.000791282  0.03863627 NaN    NaN
#> 
#> $B.mean
#> [1] 1000.0000000    0.5137792
#> 
#> $ustats
#>     ustat       sdu 
#> 0.4993414        NA 
#> 
#> attr(,"class")
#> [1] "bcaboot"

^{Created on 2021-08-26 by the reprex package (v2.0.1)}

[t-pollington]

I believe it is a problem with the size of Y sent to qbca2() as I was able to get non-NAs by making an adjustment as shown in myway2 below.

data(diabetes, package = "bcaboot")
Xy <- cbind(diabetes$x, diabetes$y)
rfun <- function(Xy) {
  y <- Xy[, 11]
  X <- Xy[, 1:10]
  summary(lm(y~X) )$adj.r.squared
}

set.seed(1234)
inds = matrix(NA, nrow = 1000, ncol = nrow(Xy))
inds = t(apply(inds, MARGIN = 1, function(inds) sample(1:nrow(Xy), replace = TRUE))) # get bootstrap indices to feed into estimator rfun()
tt = vector(mode = "numeric", length = 1000)
tt = apply(inds, MARGIN = 1, function(inds) rfun(Xy[inds,])) # get \hat{theta}*
Y = t(apply(inds, MARGIN = 1, tabulate, nbins = nrow(Xy)))

Blist = list(Y = Y, 
             tt = tt, 
             t0 = rfun(Xy)) # original estimate \hat{theta}

# bcajack2 pasted from https://github.com/bnaras/bcaboot/blob/master/R/bcajack2.R
bcaboot.return <- function(x) {
  class(x) <- "bcaboot"
  x
}
bcajack2 <- function(x, B, func, ..., m = nrow(x), mr, pct = 0.333, K = 2, J = 12,
                     alpha = c(0.025, 0.05, 0.1, 0.16),
                     verbose = TRUE) {
  
  call <- match.call()
  
  ## Save rng state
  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
    stats::runif(1)
  seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
  
  qbca2 <- function(Y, tt, t0, alpha, pct) {
    m <- ncol(Y)
    B <- nrow(Y)
    o1 <- rep(1, m)
    D <- rep(0, B)
    
    for (i in seq_len(B)) {
      Yi <- Y[i, ]
      d <- 2 * Yi * log(Yi) - 2 * (Yi - 1)
      d[Yi == 0] <- 2
      D[i] <- sum(d)
    }
    Qd <- stats::quantile(D, pct)
    ip <- seq_len(B)[D <= Qd]
    ty. <- as.vector(m * stats::lm(tt[ip] ~ Y[ip, ] - 1)$coef)
    ty. <- ty. - mean(ty.)
    a <- (1/6) * sum(ty.^3)/sum(ty.^2)^1.5
    ## if (sw == 3)
    ##     return(ty.)
    s <- mean(tt)
    B.mean <- c(B, s)
    
    zalpha <- stats::qnorm(alpha)
    nal <- length(alpha)
    ustat <- 2 * t0 - s
    ##s. <- m * .v(stats::cov(tt, Y))
    s. <- m * as.vector(stats::cov(tt, Y))
    u. <- 2 * ty. - s.
    sdu <- sum(u.^2)^0.5/m
    ustats <- c(ustat, sdu)
    names(ustats) <- c("ustat", "sdu")
    
    sdboot <- stats::sd(tt)
    sdjack <- sqrt(sum(ty.^2))/(m - 1)
    z0 <- stats::qnorm(sum(tt < t0)/B)
    
    iles <- stats::pnorm(z0 + (z0 + zalpha)/(1 - a * (z0 + zalpha)))
    ooo <- trunc(iles * B)
    ooo <- pmin(pmax(ooo, 1), B)
    lims <- sort(tt)[ooo]
    standard <- t0 + sdboot * stats::qnorm(alpha)
    ##lims <- round(cbind(lims, standard), rou)
    lims <- cbind(lims, standard)
    dimnames(lims) <- list(alpha, c("bca", "std"))
    stats <- c(t0, sdboot, z0, a, sdjack)
    names(stats) <- c("theta", "sdboot", "z0", "a", "sdjack")
    vl <- list(lims = lims, stats = stats, B.mean = B.mean, ustats = ustats)
    return(vl)
  }
  
  alpha <- alpha[alpha < 0.5]
  alpha <- c(alpha, 0.5, rev(1 - alpha))
  
  if (is.list(B)) {
    Y <- B$Y
    tt <- B$tt
    t0 <- B$t0
    B <- length(tt)
    print(paste("Y",str(Y)))
    print(paste("tt",str(tt)))
    print(paste("t0",str(t0)))
    print(paste("B",str(B)))
    print(paste("alpha",str(alpha)))
    print(paste("pct",str(pct)))
    ##vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct, rou = rou)
    vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct)
  } else {
    if (is.vector(x))
      x <- as.matrix(x)
    n <- nrow(x)
    tt <- rep(0, B)
    t0 <- func(x, ...)
    
    if (m == n) {
      ii <- sample(x = seq_len(n), size = n * B, replace = TRUE)
      ii <- matrix(ii, B)
      Y <- matrix(0, B, n)
      if (verbose) pb <- utils::txtProgressBar(min = 0, max = B, style = 3)
      for (k in 1:B) {
        ik <- ii[k, ]
        tt[k] <- func(x[ik, ], ...)
        Y[k, ] <- table(c(ik, 1:n)) - 1
        if (verbose) utils::setTxtProgressBar(pb, k)
      }
      if (verbose) close(pb)
      ##vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct, rou = rou)
      vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct)
    }
    
    if (m < n) {
      r <- n%%m
      Imat <- matrix(sample(x = seq_len(n), size = n - r), m)
      Iout <- setdiff(1:n, Imat)
      ii <- sample(x = seq_len(m), size = m * B, replace = TRUE)
      ii <- matrix(ii, B)
      Y <- matrix(0, B, m)
      if (verbose) pb <- utils::txtProgressBar(min = 0, max = B, style = 3)
      for (k in 1:B) {
        ik <- ii[k, ]
        Ik <- c(t(Imat[ik, ]))
        Ik <- c(Ik, Iout)
        tt[k] <- func(x[Ik, ], ...)
        Y[k, ] <- table(c(ik, 1:m)) - 1
        if (verbose) utils::setTxtProgressBar(pb, k)
      }
      if (verbose) close(pb)
      ##vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct, rou = rou)
      print(paste("Y",str(Y)))
      print(paste("tt",str(tt)))
      print(paste("t0",str(t0)))
      print(paste("B",str(B)))
      print(paste("alpha",str(alpha)))
      print(paste("pct",str(pct)))
      vl0 <- qbca2(Y, tt, t0, alpha = alpha, pct = pct)
    }
  }
  vl0$call <- call
  vl0$seed <- seed
  
  if (K == 0)
    bcaboot.return(vl0)
  
  nal <- length(alpha)
  Pct <- rep(0, nal)
  ##for (i in 1:nal) Pct[i] <- round(sum(tt <= vl0$lims[i, 1])/B, rou)
  for (i in 1:nal) Pct[i] <- sum(tt <= vl0$lims[i, 1])/B
  Stand <- vl0$stats[1] + vl0$stats[2] * stats::qnorm(alpha)
  Limsd <- matrix(0, length(alpha), K)
  Statsd <- matrix(0, 5, K)
  
  Limbcsd <- matrix(0, length(alpha), K)
  Statsd <- matrix(0, 5, K)
  for (k in 1:K) {
    rr <- B%%J
    II <- sample(x = B, size = B - rr)
    II <- matrix(II, ncol = J)
    limbc <- limst <- matrix(0, length(alpha), J)
    stats <- matrix(0, 5, J)
    
    for (j in 1:J) {
      iij <- c(II[, -j])
      Yj <- Y[iij, ]
      ttj <- tt[iij]
      ##vlj <- qbca2(Yj, ttj, t0, alpha, pct, rou)
      vlj <- qbca2(Yj, ttj, t0, alpha, pct)
      limbc[, j] <- vlj$lims[, 1]
      limst[, j] <- vlj$lims[, 2]
      stats[, j] <- vlj$stats
    }
    
    ## if (sw == 4)
    ##     return(list(limbc = limbc, limst = limst, stats = stats))
    Limbcsd[, k] <- apply(limbc, 1, sd) * (J - 1)/sqrt(J)
    Statsd[, k] <- apply(stats, 1, sd) * (J - 1)/sqrt(J)
    ## if (verbose)
    ##     cat("{", k, "}", sep = "")
    ## if (sw == 6)
    ##     return(list(Limbcsd = Limbcsd, Statsd = Statsd))
  }
  limsd <- rowMeans(Limbcsd, 1)
  statsd <- rowMeans(Statsd, 1)
  ##limits <- round(cbind(vl0$lims[, 1], limsd, vl0$lims[, 2], Pct), rou)
  limits <- cbind(vl0$lims[, 1], limsd, vl0$lims[, 2], Pct)
  dimnames(limits) <- list(alpha, c("bca", "jacksd", "std", "pct"))
  ##stats <- round(rbind(vl0$stats, statsd), rou)
  stats <- rbind(vl0$stats, statsd)
  ##ustats <- round(vl0$ustats, rou)
  ustats <- vl0$ustats
  ##B.mean <- c(B, round(mean(tt), rou))
  B.mean <- c(B, mean(tt))
  dimnames(stats) <- list(c("est", "jsd"), c("theta", "sdboot", "z0", "a", "sdjack"))
  vll <- list(call = call, lims = limits, stats = stats, B.mean = B.mean, ustats = ustats,
              seed = seed)
  ## if (sw == 5)
  ##     vll$tt <- tt
  bcaboot.return(vll)
}

myway = bcajack2(B = Blist, m = 40, verbose = FALSE)
#>  int [1:1000, 1:442] 0 1 0 0 0 1 2 0 1 2 ...
#> [1] "Y "
#>  num [1:1000] 0.559 0.444 0.533 0.509 0.527 ...
#> [1] "tt "
#>  num 0.507
#> [1] "t0 "
#>  int 1000
#> [1] "B "
#>  num [1:9] 0.025 0.05 0.1 0.16 0.5 0.84 0.9 0.95 0.975
#> [1] "alpha "
#>  num 0.333
#> [1] "pct "
myway$lims[,1] # Notice different size of Y sent to qbca2() vs theirway.  
#> 0.025  0.05   0.1  0.16   0.5  0.84   0.9  0.95 0.975 
#>    NA    NA    NA    NA    NA    NA    NA    NA    NA
theirway = bcajack2(x = Xy, B = 1000, func = rfun, m = 40, verbose = FALSE)
#>  num [1:1000, 1:40] 2 2 0 0 0 2 0 1 1 0 ...
#> [1] "Y "
#>  num [1:1000] 0.543 0.605 0.551 0.525 0.462 ...
#> [1] "tt "
#>  num 0.507
#> [1] "t0 "
#>  num 1000
#> [1] "B "
#>  num [1:9] 0.025 0.05 0.1 0.16 0.5 0.84 0.9 0.95 0.975
#> [1] "alpha "
#>  num 0.333
#> [1] "pct "
theirway$lims[,1]
#>     0.025      0.05       0.1      0.16       0.5      0.84       0.9      0.95 
#> 0.4119962 0.4243488 0.4469171 0.4568380 0.4952950 0.5278410 0.5380414 0.5497895 
#>     0.975 
#> 0.5609942

myway2 = bcajack2(B = list(Y = Y[,1:40], tt = tt, t0 = rfun(Xy)), m = 40, 
                  verbose = FALSE)
#>  int [1:1000, 1:40] 0 1 0 0 0 1 2 0 1 2 ...
#> [1] "Y "
#>  num [1:1000] 0.559 0.444 0.533 0.509 0.527 ...
#> [1] "tt "
#>  num 0.507
#> [1] "t0 "
#>  int 1000
#> [1] "B "
#>  num [1:9] 0.025 0.05 0.1 0.16 0.5 0.84 0.9 0.95 0.975
#> [1] "alpha "
#>  num 0.333
#> [1] "pct "
myway2$lims[,1]
#>     0.025      0.05       0.1      0.16       0.5      0.84       0.9      0.95 
#> 0.4355705 0.4436322 0.4566175 0.4664901 0.5002812 0.5326244 0.5433087 0.5533786 
#>     0.975 
#> 0.5619291

^{Created on 2021-08-26 by the reprex package (v2.0.1)}

[t-pollington]

I'm more confident of the bug. It is not due to line 98 producing -Infs. It is caused on line 104 by stats::lm(y~x) requiring that there are more/equal observations than the number of dimensions of the design matrix. The lm()$coeffs term will contain NAs when length(ip)(no. of obs) < ncol(Y) (no. of dimensions p, where p I mean the n x p design matrix (x) as defined internally in stats::lm.fit(x,...)).

Furthermore this will occur when:
length(ip) < ncol(Y)
or equivalently B x pct (approximately) <= ncol(Y) = n
I say 'approximately' because it depends on the random sample of bootstrapped indices in each Y[i,] which affects the amount of D<Qd.

Herein lies the problem:

1. For the default pct=0.33 and B=2500 then if n is approximately equal or bigger than 825 then lm() will be unable to produce an NA-free vector. This limits bcajack2() whose advantage is not having to do n replications to estimate a.
2. In my case where I have n=25000, even if I had done more than 75K bootstrap replications, lm() would run out of typical available RAM since the design matrix x would need to be 25000x25000.

Thus bcajack() seems to be the best bet for dataset sizes of 800 and above. I plan to submit a short pull request that errors when lm() silently produces this NA-containing vector.

[bnaras]

Hi @t-pollington , thank you for your detailed reports and the PR. I will have a chance to look at it in detail in the next few days.

[t-pollington]

Glad it's helpful and thanks for sharing this package. I still don't understand why BCa CIs aren't as ubiquitous in data analysis as AIC.
Thank you very much for a well written article which links seamlessly with this code.