neg_bin_pi()
calculates bootstrap calibrated prediction intervals for
negative-binomial data.
neg_bin_pi(
histdat,
newdat = NULL,
newoffset = 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"
)
neg_bin_pi()
returns an object of class c("predint", "negativeBinomialPI")
with prediction intervals or limits in the first entry ($prediction
).
a data.frame
with two columns. The first has to contain
the historical observations. The second has to contain the number of experimental
units per study (offsets).
data.frame
with two columns. The first has to contain
the future observations. The second has to contain the number of experimental
units per study (offsets).
vector with future number of experimental units per historical study.
either "both", "upper" or "lower".
alternative
specifies if a prediction interval or
an upper or a lower prediction limit should be computed
defines the level of confidence (\(1-\alpha\))
number of bootstraps
lower start value for bisection
upper start value for bisection
tolerance for the coverage probability in the bisection
if TRUE
: Plot for visualization of the bisection process
maximal number of bisection steps
either "MS22" or "MS22mod" (see 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 = n^*_m \hat{\lambda} \pm q \sqrt{n^*_m \frac{\hat{\lambda} + \hat{\kappa} \bar{n} \hat{\lambda}}{\bar{n} H} + (n^*_m \hat{\lambda} + \hat{\kappa} n^{*2}_m \hat{\lambda}^2) }$$
with \(n^*_m\) as the number of experimental units in the future clusters,
\(\hat{\lambda}\) as the estimate for the Poisson mean obtained from the
historical data, \(\hat{\kappa}\) as the estimate for the dispersion parameter,
\(n_h\) as the number of experimental units per historical cluster and
\(\bar{n}=\sum_h^{n_h} n_h / H\).
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] = \Big[n^*_m \hat{\lambda} - q^{calib}_l \sqrt{n^*_m \frac{\hat{\lambda} + \hat{\kappa} \bar{n} \hat{\lambda}}{\bar{n} H} + (n^*_m \hat{\lambda} + \hat{\kappa} n^{*2}_m \hat{\lambda}^2)}, \quad n^*_m \hat{\lambda} + q^{calib}_u \sqrt{n^*_m \frac{\hat{\lambda} + \hat{\kappa} \bar{n} \hat{\lambda}}{\bar{n} H} + (n^*_m \hat{\lambda} + \hat{\kappa} n^{*2}_m \hat{\lambda}^2) } \Big]$$
Please note, that this modification does not affect the calibration procedure, if only prediction limits are of interest.
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")
# HCD from the Ames test
ames_HCD
# Prediction interval for one future number of revertant colonies
# obtained in three petridishes
pred_int <- neg_bin_pi(histdat=ames_HCD, newoffset=3, nboot=100)
summary(pred_int)
# Please note that nboot was set to 100 in order to decrease computing time
# of the example. For a valid analysis set nboot=10000.
Run the code above in your browser using DataLab