beta_bin_pi function

Prediction intervals for beta-binomial data

Prediction intervals for beta-binomial data

beta_bin_pi() calculates bootstrap calibrated prediction intervals for beta-binomial data

beta_bin_pi(
  histdat,
  newdat = NULL,
  newsize = NULL,
  alternative = "both",
  alpha = 0.05,
  nboot = 10000,
  delta_min = 0.01,
  delta_max = 10,
  tolerance = 0.001,
  traceplot = TRUE,
  n_bisec = 30,
  algorithm = "MS22mod"
)

Arguments

  • histdat: a data.frame with two columns (number of successes and number of failures) containing the historical data

  • newdat: a data.frame with two columns (number of successes and number of failures) containing the future data

  • newsize: a vector containing the future cluster sizes

  • alternative: either "both", "upper" or "lower". alternative

    specifies if a prediction interval or an upper or a lower prediction limit should be computed

  • alpha: defines the level of confidence (1-alpha)

  • nboot: number of bootstraps

  • delta_min: lower start value for bisection

  • delta_max: upper start value for bisection

  • tolerance: tolerance for the coverage probability in the bisection

  • traceplot: if TRUE: Plot for visualization of the bisection process

  • n_bisec: maximal number of bisection steps

  • algorithm: either "MS22" or "MS22mod" (see details)

Returns

beta_bin_pi returns an object of class c("predint", "betaBinomialPI")

with prediction intervals or limits in the first entry ($prediction).

Details

This function returns bootstrap-calibrated prediction intervals as well as lower or upper prediction limits.

If algorithm is set to "MS22", both limits of the prediction interval are calibrated simultaneously using the algorithm described in Menssen and Schaarschmidt (2022), section 3.2.4. The calibrated prediction interval is given as

[l,u]m=nm∗π^±qcalibse^(Ym−ym∗) [l,u]_m = n^*_m \hat{\pi} \pm q^{calib} \hat{se}(Y_m - y^*_m)

where

se^(Ym−ym∗)=nm∗π^(1−π^)[1+(nm∗−1)ρ^]+[nm∗2π^(1−π^)∑hnh+∑hnh−1∑hnhnm∗2π^(1−π^)ρ^] \hat{se}(Y_m - y^*_m) = \sqrt{n^*_m \hat{\pi} (1- \hat{\pi}) [1 + (n^*_m -1) \hat{\rho}] +[\frac{n^{*2}_m \hat{\pi} (1- \hat{\pi})}{\sum_h n_h} + \frac{\sum_h n_h -1}{\sum_h n_h} n^{*2}_m \hat{\pi} (1- \hat{\pi}) \hat{\rho}]}

with nm∗n^*_m as the number of experimental units in the future clusters, π^\hat{\pi} as the estimate for the binomial proportion obtained from the historical data, qcalibq^{calib} as the bootstrap-calibrated coefficient, ρ^\hat{\rho} as the estimate for the intra class correlation (Lui et al. 2000) and nhn_h as the number of experimental units per historical cluster.

If algorithm is set to "MS22mod", both limits of the prediction interval are calibrated independently from each other. The resulting prediction interval is given by

[l,u]m=[nm∗π^−qlcalibse^(Ym−ym∗),nm∗π^+qucalibse^(Ym−ym∗)] [l,u]_m = \big[n^*_m \hat{\pi} - q^{calib}_l \hat{se}(Y_m - y^*_m), \quad n^*_m \hat{\pi} + q^{calib}_u \hat{se}(Y_m - y^*_m) \big]

Please note, that this modification does not affect the calibration procedure, if only prediction limits are of interest.

Examples

# Prediction interval
pred_int <- beta_bin_pi(histdat=mortality_HCD, newsize=40, nboot=100)
summary(pred_int)

# Upper prediction bound
pred_u <- beta_bin_pi(histdat=mortality_HCD, newsize=40, alternative="upper", nboot=100)
summary(pred_u)

# Please note that nboot was set to 100 in order to decrease computing time
# of the example. For a valid analysis set nboot=10000.

References

Lui et al. (2000): Confidence intervals for the risk ratio under cluster sampling based on the beta-binomial model. Statistics in Medicine.

tools:::Rd_expr_doi("10.1002/1097-0258(20001115)19:21<2933::AID-SIM591>3.0.CO;2-Q")

Menssen and Schaarschmidt (2022): Prediction intervals for all of M future observations based on linear random effects models. Statistica Neerlandica. tools:::Rd_expr_doi("10.1111/stan.12260")

predint packageRead PDF manual
  • Maintainer: Max Menssen
  • License: GPL (>= 2)
  • Last published: 2024-03-04

Useful links