boot.pval throws an error for BCa confidence intervals

[astaroph]

Just a small note to draw your attention to an error I got when calculating p-values for BCa confidence intervals along with my temporary fix (though I still don't understand why it works). Hopefully you might!

When using the original boot.pval function I get the following error when type='bca' but not when type is left to defaults:

Error in if (any(ints)) out[inds[ints]] <- tstar[k[inds[ints]]] :
missing value where TRUE/FALSE needed
6.
norm.inter(t, adj.alpha)
5.
bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o,
h = h, hdot = hdot, hinv = hinv, ...)
4.
boot::boot.ci(boot_res, conf = 1 - alpha_seq, type = type, ...)
3.
withCallingHandlers(expr, warning = function(w) if (inherits(w,
classes)) tryInvokeRestart("muffleWarning"))
2.
suppressWarnings(boot::boot.ci(boot_res, conf = 1 - alpha_seq,
type = type, ...))
1.
boot.pval_fixed(boot_out, type = "bca", theta_null = 0.5)

I apparently fixed it by changing a single value in the following line:
original (doesn't work with type='bca'):
alpha_seq <- seq(1e-16, 1-1e-16, pval_precision)

modified (does work with type='bca')
alpha_seq <- seq(1e-15, 1-1e-15, pval_precision)

Strange!

[mthulin]

Thanks for posting this! I'll look into it as soon as I can.

[astaroph]

Sure thing!

I also noticed something else:
When running this again on another dataset
alpha_seq <- seq(1e-15, 1-1e-15, pval_precision)
also failed, but
alpha_seq <- seq(pval_precision, 1-pval_precision, pval_precision)
seems to work every time (so far), although this now has the downside of limiting the minimum possible p-value by the number of bootstraps...
Anyhow, just wanted to give you a heads up, but thanks for putting out there a truly useful solution for generating p-values from bootstrapping results!

[arielhasidim]

+1

cluster_boot_out <- boot(pair_ids, cluster_boot_fun, R = 9999)
boot.pval(cluster_boot_out, 
          theta_null = 1,
          type = "bca")

# Error in if (any(ints)) out[inds[ints]] <- tstar[k[inds[ints]]] :
# missing value where TRUE/FALSE needed
# 6.
# norm.inter(t, adj.alpha)
# 5.
# bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o,
# h = h, hdot = hdot, hinv = hinv, ...)
# 4.
# boot::boot.ci(boot_res, conf = 1 - alpha_seq, type = type, ...)
# 3.
# withCallingHandlers(expr, warning = function(w) if (inherits(w,
# classes)) tryInvokeRestart("muffleWarning"))
# 2.
# suppressWarnings(boot::boot.ci(boot_res, conf = 1 - alpha_seq,
# type = type, ...))
# 1.
# boot.pval(cluster_boot_out, theta_null = 1, type = "bca")

[Sciurus365]

Hi, I had the same problem and I figured out why it happened. In boot.pval(), the alpha sequence is set as alpha_seq <- seq(1e-16, 1-1e-16, pval_precision), and 1-alpha_seq is eventually passed to boot:::bca.ci() as the confidence level. But 1e-16 is too small that although R thinks 1e-16!=0, it thinks (1+1-1e-16)/2==1, and qnorm((1+1-1e-16)/2) is Inf (a step in boot:::bca.ci()). (So if a value's difference to 1 is smaller than 5e-17, R thinks it is exactly 1. You can test it by checking 1-5e-17 == 1.) This Inf later leads to the missing value reported above.
Making 1e-16 a little bit larger will solve this problem, as suggested by @astaroph. I'll submit a pull request if that helps. @mthulin
(I didn't encounter any cases where changing it to 1e-15 still leads to problems, but I'm also happy to look into it. @astaroph In case you are still interested in this issue and have a reproducible example, just let me know!)

[astaroph]

@Sciurus365 It's been a little bit since I revisited this issue and my datasets have evolved a bit since then, so I couldn't find the specific example where setting it to 1e-15 failed (though I know it did). However, I think I may have traced the issue to the pval_precision parameter. Observe:
Assuming you are defining the alpha_seq the same way and all you change is the pval_precision parameter (which will vary with the number of boots as it is defined as pval_precision = 1/boot_res$R if left blank then:
alpha_seq=seq(1e-15, 1-1e-15, 1e-05) will yield a proper sequence of 100000 elements ending with 0.99999
However alpha_seq=seq(1e-15, 1-1e-15, 1e-04) will yield a sequence of 10,001 elements the last of which is a NA, yet
alpha_seq=seq(1e-15, 1-1e-15, 1e-06) will yield a proper sequence of 1,000,000 elements ending with 0.999999.

This suggests that the step size of the sequence set by pval_precision (and which depends on the number of boot supplied by the user) interacts with the specific minimum value set in the alpha sequence in a way that does not completely depend on how large this minimum value is. Pragmatically, I could not reproduce this error trying a range of different precision values of different orders of magnitude when I used 1e-14 as opposed to 1e-15, so perhaps this is just an edge case that only happens when the step size and sequence boundary sizes are just small enough to trigger a NA at the end of the sequence? I will try setting mine to 1e-14 and see if anything ever triggers it and report back if it does.

[arielhasidim]

+1 to @astaroph. I also got the error when setting it to 1e-15.

This is my workaround (inside, there is a link to the patchwork code):
https://stats.stackexchange.com/questions/635200/can-i-remove-a-single-bootstrap-sample-before-calculating-ci

[Sciurus365]

@Sciurus365 It's been a little bit since I revisited this issue and my datasets have evolved a bit since then, so I couldn't find the specific example where setting it to 1e-15 failed (though I know it did). However, I think I may have traced the issue to the pval_precision parameter. Observe: Assuming you are defining the alpha_seq the same way and all you change is the pval_precision parameter (which will vary with the number of boots as it is defined as pval_precision = 1/boot_res$R if left blank then: alpha_seq=seq(1e-15, 1-1e-15, 1e-05) will yield a proper sequence of 100000 elements ending with 0.99999 However alpha_seq=seq(1e-15, 1-1e-15, 1e-04) will yield a sequence of 10,001 elements the last of which is a NA, yet alpha_seq=seq(1e-15, 1-1e-15, 1e-06) will yield a proper sequence of 1,000,000 elements ending with 0.999999.

This suggests that the step size of the sequence set by pval_precision (and which depends on the number of boot supplied by the user) interacts with the specific minimum value set in the alpha sequence in a way that does not completely depend on how large this minimum value is. Pragmatically, I could not reproduce this error trying a range of different precision values of different orders of magnitude when I used 1e-14 as opposed to 1e-15, so perhaps this is just an edge case that only happens when the step size and sequence boundary sizes are just small enough to trigger a NA at the end of the sequence? I will try setting mine to 1e-14 and see if anything ever triggers it and report back if it does.

Hey, thanks for your reply! Do you mean the vector from seq(1e-15, 1-1e-15, 1e-04) yields an NA at the end, or it leads to an NA later in the code? Because I just tried seq(1e-15, 1-1e-15, 1e-04) and the last element equals 1-1e-15 and is not NA.
(If you get an NA here directly, which R version are you using?)