| Title: | Inference with Extreme Value Inflated Count Data |
|---|---|
| Description: | Allows users to model and draw inferences from extreme value inflated count data, and to evaluate these models and compare them to their non extreme-value inflated counterparts. The package is built to be compatible with standard model-presentation tools such as the 'broom', 'generics', and 'modelsummary' packages. |
| Authors: | David Randahl [cre, aut], Johan Vegelius [aut] |
| Maintainer: | David Randahl <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.0 |
| Built: | 2026-09-24 19:41:04 UTC |
| Source: | https://github.com/doktorandahl/evinf |
Fits n additional bootstrap replicates with the same block structure
and appends them, so you don't have to refit everything when the original
number turns out to be too few.
add_bootstraps( object, n, boot_seed = NULL, multicore = NULL, ncores = NULL, verbose = FALSE )add_bootstraps( object, n, boot_seed = NULL, multicore = NULL, ncores = NULL, verbose = FALSE )
object |
A fitted |
n |
Number of additional bootstrap replicates. |
boot_seed |
RNG seed for this batch. Each batch needs its own seed:
the bootstrap stream is fully determined by the seed, so reusing one that
the model (or an earlier |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
verbose |
Show a progress bar. |
The model with n more replicates in object$bootstraps
(names continue bootstrap_<k>).
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) model <- add_bootstraps(model, 5)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) model <- add_bootstraps(model, 5)
Augment data with fitted values, residuals and state classification
## S3 method for class 'evinb' augment(x, data = NULL, newdata = NULL, seed = NULL, ...)## S3 method for class 'evinb' augment(x, data = NULL, newdata = NULL, seed = NULL, ...)
x |
A fitted |
data, newdata
|
Data to augment. |
seed |
Optional seed for the randomized quantile residual
( |
... |
Unused. |
See augment.evzinb; evinb has no zero state,
so .prob_zero/.post_zero are not produced.
data(genevzinb2) model <- evinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) augment(model)data(genevzinb2) model <- evinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) augment(model)
Augment data with fitted values, residuals and state classification
## S3 method for class 'evzinb' augment(x, data = NULL, newdata = NULL, seed = NULL, ...)## S3 method for class 'evzinb' augment(x, data = NULL, newdata = NULL, seed = NULL, ...)
x |
A fitted |
data, newdata
|
Data to augment. |
seed |
Optional seed for the randomized quantile residual
( |
... |
Unused. |
data/newdata (or the estimation data) with columns
appended: .fitted (the predict(type = "harmonic") point
prediction), .prob_zero/.prob_count/.prob_evi (prior
state probabilities – augment.evinb() has no zero state, so
.prob_zero is dropped there) and .state (the MAP prior
state, from classify_states()). When the response column is
present, three more columns are added: .resid (a seeded
randomized quantile residual, residuals(type = "quantile")),
.post_zero/.post_count/.post_evi (posterior state
probabilities) and .post_state (the MAP posterior state).
data(genevzinb2) model <- evzinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) augment(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) augment(model)
Returns the log-likelihood evaluated at each candidate value of
from the final EM iteration. Use it to judge whether is well
identified and whether the candidate range (c.lim) is wide enough.
c_profile(object)c_profile(object)
object |
A fitted |
A tibble with columns c and loglik.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) c_profile(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) c_profile(model)
A single report gathering the convergence and health checks otherwise
spread across print(), glance() and failed_bootstraps():
EM and -profile convergence, whether landed on a
candidate-grid boundary (full sample and across bootstrap replicates), the
smallest fitted Pareto shape against alpha_pl_floor, failed and
degenerate bootstrap replicates, start agreement (see n_starts in
evinf_control), and, for the block bootstrap schemes, the
out-of-bag row fraction.
check_evinf(object)check_evinf(object)
object |
A fitted |
A classed tibble ("evinf_check") with columns check,
status ("ok", "note" or "warning") and
detail (a one-line description, with a suggestion for anything
not "ok").
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) check_evinf(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) check_evinf(model)
For each observation, returns the prior state probabilities (from the inflation components) and, for the estimation data, the posterior state probabilities (responsibilities), together with the maximum-a-posteriori state under each.
classify_states( object, rule = c("map", "threshold"), threshold = 0.5, newdata = NULL )classify_states( object, rule = c("map", "threshold"), threshold = 0.5, newdata = NULL )
object |
A fitted |
rule |
|
threshold |
Threshold for |
newdata |
Optional data. Prior probabilities are always available;
posterior probabilities need the response column and are |
A tibble with prior_zero / prior_count / prior_evi,
posterior_*, map_prior, map_posterior (factors with
levels zero, count, evi) and y.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) classify_states(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) classify_states(model)
Coefficients of an evzinb / evinb model
## S3 method for class 'evzinb' coef(object, component = "all", ...) ## S3 method for class 'evinb' coef(object, component = "all", ...)## S3 method for class 'evzinb' coef(object, component = "all", ...) ## S3 method for class 'evinb' coef(object, component = "all", ...)
object |
A fitted |
component |
One of |
... |
Unused. |
For "all", a named numeric vector <component>_<term>,
alpha_nb, c_ev; for a single component the plain named vector.
Bootstrap coefficient extractor
coefficient_extractor(object, ...) ## S3 method for class 'evzinb' coefficient_extractor( object, component = c("all", "count", "zero", "evi", "pareto"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evinb' coefficient_extractor( object, component = c("all", "count", "evi", "pareto"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'zinbboot' coefficient_extractor(object, component = c("all", "count", "zero"), ...) ## S3 method for class 'nbboot' coefficient_extractor(object, ...)coefficient_extractor(object, ...) ## S3 method for class 'evzinb' coefficient_extractor( object, component = c("all", "count", "zero", "evi", "pareto"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evinb' coefficient_extractor( object, component = c("all", "count", "evi", "pareto"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'zinbboot' coefficient_extractor(object, component = c("all", "count", "zero"), ...) ## S3 method for class 'nbboot' coefficient_extractor(object, ...)
object |
a fitted model with bootstraps of class evzinb, evinb, nbboot, or zinbboot |
... |
Arguments passed to methods, in particular |
component |
Which component should be extracted |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate
(default |
A tibble with coefficient values, one row per bootstrap and component
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) coefficient_extractor(model, component = 'all') data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) coefficient_extractor(zinb_comp$zinb) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) coefficient_extractor(zinb_comp$nb)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) coefficient_extractor(model, component = 'all') data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) coefficient_extractor(zinb_comp$zinb) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) coefficient_extractor(zinb_comp$nb)
For each compared model (nb, zinb, and any winsorized /
razorized variants) and each metric, computes the paired bootstrap difference
evinf - compared (so a negative median favours the extreme-value
model, matching Table B3 of the ISQ appendix), the proportion of bootstraps
in which the evinf model is better, and the number of bootstrap pairs where
both fits succeeded.
compare_fit( comp, metrics = c("aic", "bic", "rmse", "rmsle"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evzinbcomp' plot( x, metrics = c("aic", "bic", "rmse", "rmsle"), exclude_degenerate = TRUE, ... )compare_fit( comp, metrics = c("aic", "bic", "rmse", "rmsle"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evzinbcomp' plot( x, metrics = c("aic", "bic", "rmse", "rmsle"), exclude_degenerate = TRUE, ... )
comp |
An |
metrics |
Metrics to show. |
exclude_degenerate |
Drop bootstrap replicates of the evinf model
flagged degenerate (default |
... |
Unused. |
x |
An |
The aic / bic rows are NA for a *_razor or
*_winsor slot (round8 0.5): a razorised fit is estimated on fewer
observations than the evinf model, and a winsorised fit is estimated on a
different outcome, so their AIC/BIC are not on the same scale as the evinf
model's and a difference between them is not meaningful. The RMSE / RMSLE
rows for those slots remain comparable, because out-of-bag error is always
computed against the raw (un-winsorised, un-razorised) outcome.
A tibble of class evinf_compare_fit with columns model,
metric, median_difference, prop_evinf_better,
n_pairs. median_difference / prop_evinf_better are
NA for the aic / bic rows of a *_razor or
*_winsor slot (see Details); n_pairs is left as computed on
the returned object (print() blanks it for those rows instead,
since a pair count next to an NA metric reads like a bug).
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 10) compare_fit(compare_models(model))data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 10) compare_fit(compare_models(model))
Function to compare evzinb or evinb models with zinb and nb models
compare_models( object, nb_comparison = TRUE, zinb_comparison = TRUE, poisson_comparison = identical((object$family %||% evinf_family())$count, "poisson"), zip_comparison = identical((object$family %||% evinf_family())$count, "poisson"), winsorize = FALSE, razorize = FALSE, cutoff_value = 10, init_theta = NULL, multicore = NULL, ncores = NULL )compare_models( object, nb_comparison = TRUE, zinb_comparison = TRUE, poisson_comparison = identical((object$family %||% evinf_family())$count, "poisson"), zip_comparison = identical((object$family %||% evinf_family())$count, "poisson"), winsorize = FALSE, razorize = FALSE, cutoff_value = 10, init_theta = NULL, multicore = NULL, ncores = NULL )
object |
A fitted evzinb or evinb model object |
nb_comparison |
Should comparison be made with a negative binomial model? |
zinb_comparison |
Should comparisons be made with the zinb model? Not
available for |
poisson_comparison |
Should comparison be made with a Poisson model?
Defaults to |
zip_comparison |
Should comparisons be made with a zero-inflated
Poisson (ZIP) model? Same default as |
winsorize |
Should winsorizing be done in the comparisons? |
razorize |
Should razorizing (trimming) be done in the comparisons? |
cutoff_value |
Integer: Which observation should be used as a basis for winsorizing/razorising. E.g. 10 means that everything larger than the 10th observation will be winsorized/razorised |
init_theta |
Optional initial value for theta in the NB specification |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
An object of class evzinbcomp: a list whose first element
model is the original evzinb/evinb model (also available as
evzinb for backwards compatibility), followed by the compared
nb / zinb models (and their winsorized/razorized variants when
requested).
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) compare_models(model) data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log + osvAllLagDum, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE) cmp <- compare_models(hks_mod) compare_fit(cmp)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) compare_models(model) data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log + osvAllLagDum, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE) cmp <- compare_models(hks_mod) compare_fit(cmp)
Confidence intervals for an evzinb / evinb model
## S3 method for class 'evzinb' confint( object, parm, level = 0.95, type = c("percentile", "approx"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evinb' confint( object, parm, level = 0.95, type = c("percentile", "approx"), exclude_degenerate = TRUE, ... )## S3 method for class 'evzinb' confint( object, parm, level = 0.95, type = c("percentile", "approx"), exclude_degenerate = TRUE, ... ) ## S3 method for class 'evinb' confint( object, parm, level = 0.95, type = c("percentile", "approx"), exclude_degenerate = TRUE, ... )
object |
A fitted model with bootstraps. |
parm |
Which parameters (names as in |
level |
Confidence level. |
type |
|
exclude_degenerate |
Drop bootstrap replicates flagged degenerate
(default |
... |
Unused. |
A two-column matrix.
data(genevzinb2) m <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 25) # the Pareto shape (alpha_nb) and the threshold (c_ev) are in coef()/confint() confint(m, parm = c("count_x1", "alpha_nb", "c_ev"))data(genevzinb2) m <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 25) # the Pareto shape (alpha_nb) and the threshold (c_ev) are in coef()/confint() confint(m, parm = c("count_x1", "alpha_nb", "c_ev"))
Running an extreme value inflated negative binomial model with bootstrapping
evinb( formula_nb, formula_evi = NULL, formula_pareto = NULL, data, bootstrap = TRUE, n_bootstraps = 100, multicore = NULL, ncores = NULL, block = NULL, weights = NULL, time = NULL, bootstrap_scheme = NULL, block_length = NULL, family = evinf_family(), boot_seed = NULL, start_seed = NULL, control = evinf_control(), max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C, verbose = FALSE )evinb( formula_nb, formula_evi = NULL, formula_pareto = NULL, data, bootstrap = TRUE, n_bootstraps = 100, multicore = NULL, ncores = NULL, block = NULL, weights = NULL, time = NULL, bootstrap_scheme = NULL, block_length = NULL, family = evinf_family(), boot_seed = NULL, start_seed = NULL, control = evinf_control(), max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C, verbose = FALSE )
formula_nb |
Formula for the negative binomial (count) component of the model.
May include an |
formula_evi |
Formula for the extreme-value inflation component of the model.
If NULL taken as the same formula as nb, with any |
formula_pareto |
Formula for the pareto (extreme value) component of the
model. If NULL taken as the same formula as nb, offset stripped as above.
|
data |
Data to run the model on |
bootstrap |
Should bootstrapping be performed. Needed to obtain standard errors and p-values |
n_bootstraps |
Number of bootstraps to run. For use of bootstrapped p-values, at least 1,000 bootstraps are recommended. For approximate p-values, a lower number can be sufficient |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
block |
Optional case-identifier column for block bootstrapping, given
either as a bare column name ( |
weights |
Optional observation weights (round9 D.2), given as a bare
column name ( |
time |
Optional time index for panel/time-series bootstrap resampling
(round9 F), given as a bare column name ( |
bootstrap_scheme |
One of |
block_length |
Block length for |
family |
An |
boot_seed |
Optional bootstrap seed for reproducibility. When supplied
it is used as-is; when |
start_seed |
Optional seed for the perturbed starts when
|
control |
An |
max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C
|
Deprecated. Pass these through |
verbose |
Should progress be printed for the first run of evinb |
An object of class 'evinb'
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) data(hks) hks_mod <- evinb( osvAll ~ troopLag + policeLag + militaryobserversLag + epduration + lntpop + brv_AllLag_log + osvAllLagDum + incomp, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE )data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) data(hks) hks_mod <- evinb( osvAll ~ troopLag + policeLag + militaryobserversLag + epduration + lntpop + brv_AllLag_log + osvAllLagDum + incomp, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE )
Generates n observations from an EVZINB mixture with fixed,
hand-chosen true parameters (see Details) – formula_nb = y ~ x1 +
x2 + x3, formula_zi/formula_evi sharing the same three
covariates, formula_pareto depending on x1 alone. Meant for
benchmarking (inst/bench/bench_evinf.R) and scale testing, not as a
realistic applied example – see genevzinb2 for that.
evinf_bench_data(n = 1000, seed = NULL)evinf_bench_data(n = 1000, seed = NULL)
n |
Number of observations. |
seed |
Optional seed; the caller's |
True parameters (on each component's own scale, the same
parameterisation coef() returns):
Zero-state log-odds against the count-state baseline:
.
Evi-state log-odds against the count-state baseline:
(a deliberately small share, matching
an extreme-value-inflation process being the rare state).
Count state (negative binomial, log link):
, dispersion
.
Evi state (the discretised Pareto tail the package fits
throughout): threshold , shape (kept well above 0 for any simulated
x1, so the tail never collapses).
A tibble with columns y, x1, x2, x3.
d <- evinf_bench_data(200, seed = 1) m <- evzinb(y ~ x1 + x2 + x3, data = d, bootstrap = FALSE)d <- evinf_bench_data(200, seed = 1) m <- evzinb(y ~ x1 + x2 + x3, data = d, bootstrap = FALSE)
Bundles the ~20 tuning arguments of the EM fitting routine into a single
object, in the style of stats::glm.control() or
pscl::zeroinfl.control(). Pass the result as the control
argument of evzinb() / evinb().
evinf_control( max.diff.par = 0.01, max.no.em.steps = 500, max.no.em.steps.warmup = 5, c.lim = NULL, prune.c.range = FALSE, max.upd.par.zc.multinomial = 0.5, max.upd.par.pl.multinomial = 0.5, max.upd.par.nb = 0.5, max.upd.par.pl = 0.5, no.m.bfgs.steps.multinomial = 3, no.m.bfgs.steps.nb = 3, no.m.bfgs.steps.pl = 3, pdf.pl.type = c("approx", "exact"), eta.int = c(-1, 1), init.Beta.multinom.ZC = NULL, init.Beta.multinom.PL = NULL, init.Beta.NB = NULL, init.Beta.PL = NULL, init.Alpha.NB = 0.01, init.C = NULL, alpha_floor = 0.001, coef_limit = 50, max.c.iter = 50, alpha_pl_floor = 0.01, n_starts = 1L, start_jitter = 0.5, chunk_size = NULL )evinf_control( max.diff.par = 0.01, max.no.em.steps = 500, max.no.em.steps.warmup = 5, c.lim = NULL, prune.c.range = FALSE, max.upd.par.zc.multinomial = 0.5, max.upd.par.pl.multinomial = 0.5, max.upd.par.nb = 0.5, max.upd.par.pl = 0.5, no.m.bfgs.steps.multinomial = 3, no.m.bfgs.steps.nb = 3, no.m.bfgs.steps.pl = 3, pdf.pl.type = c("approx", "exact"), eta.int = c(-1, 1), init.Beta.multinom.ZC = NULL, init.Beta.multinom.PL = NULL, init.Beta.NB = NULL, init.Beta.PL = NULL, init.Alpha.NB = 0.01, init.C = NULL, alpha_floor = 0.001, coef_limit = 50, max.c.iter = 50, alpha_pl_floor = 0.01, n_starts = 1L, start_jitter = 0.5, chunk_size = NULL )
max.diff.par |
EM convergence tolerance: the algorithm has converged when
the maximum absolute change in the parameter estimates falls below this.
The EM log-likelihood can have close local optima; a tighter tolerance
reduces (without eliminating) the chance that two future backends run
with the same |
max.no.em.steps |
Maximum number of EM steps. |
max.no.em.steps.warmup |
Number of EM steps in each warm-up round. |
c.lim |
|
prune.c.range |
|
max.upd.par.zc.multinomial, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl
|
Maximum parameter-change step sizes for the zero-inflation, extreme-value inflation, count and Pareto components. |
no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl
|
Number of BFGS steps per M-step for the multinomial, count and Pareto blocks. |
pdf.pl.type |
Which Pareto-block derivatives the M-step uses for
|
eta.int |
Interval for the eta line search, a numeric vector of length 2. |
init.Beta.multinom.ZC, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL
|
Optional starting values for the component coefficient vectors ( |
init.Alpha.NB |
Starting value for the negative-binomial dispersion. |
init.C |
|
alpha_floor, coef_limit
|
Thresholds for flagging a bootstrap replicate
as degenerate ( |
max.c.iter |
Maximum number of outer ECME iterations (each one re-profiles
|
alpha_pl_floor |
Floor for the fitted Pareto shape |
n_starts |
Number of starting points for the full-sample fit (round10
G.1, audit §5.10). |
start_jitter |
Standard deviation of the |
chunk_size |
Number of bootstrap replicates (or starts) handed to a
worker at a time in |
When c.lim = NULL, evzinb() / evinb() choose a
range from the data: the lower bound is the 90th percentile of the positive
response values (rounded down to an observed value), the upper bound is the
third-largest unique response value; if fewer than 10 unique positive values
lie in that range the lower bound drops to the 75th percentile. The chosen
range is printed with a message and stored in object$control$c.lim.
A list of class "evinf_control".
evinf_control(max.no.em.steps = 300, c.lim = c(20, 500))evinf_control(max.no.em.steps = 300, c.lim = c(20, 500))
Specifies the distributional family of the count state (count) and
of the zero state (zero). One object, so later additions (a
truncated-count state, a geometric count) do not add more arguments to
evzinb() / evinb().
evinf_family(count = c("nbinom", "poisson"), zero = c("mixture", "hurdle"))evinf_family(count = c("nbinom", "poisson"), zero = c("mixture", "hurdle"))
count |
|
zero |
|
Sign convention against pscl (round10 0.9,
review §7): for zero = "hurdle", hurdle's
zero-hurdle component models , while evinf's zero-state
coefficients model directly – so a fitted
evzinb(family = evinf_family(zero = "hurdle")) model's zero
coefficients (Beta.multinom.ZC) are the exact negatives of the
corresponding pscl::hurdle(..., dist = "negbin") zero
coefficients, everything else about the fit being equal (verified in
test-hurdle.R). This is a clean reparameterisation because a
hurdle's zero state owns every zero outright. zero = "mixture"
(the default) has no such simple relationship to
zeroinfl's zero-inflation coefficients: a mixture's
zeros can also arise from the count state, so the two models' zero
components are not a reparameterisation of each other.
An object of class "evinf_family": a list with elements
count and zero.
evinf_family() evinf_family(count = "poisson")evinf_family() evinf_family(count = "poisson")
logLik / nobs / model.frame / formula / terms for evzinb / evinb models
## S3 method for class 'evzinb' logLik(object, ...) ## S3 method for class 'evinb' logLik(object, ...) ## S3 method for class 'evzinb' nobs(object, ...) ## S3 method for class 'evinb' nobs(object, ...) ## S3 method for class 'evzinb' model.frame(formula, ...) ## S3 method for class 'evinb' model.frame(formula, ...) ## S3 method for class 'evzinb' formula(x, component = "count", ...) ## S3 method for class 'evinb' formula(x, component = "count", ...) ## S3 method for class 'evzinb' terms(x, component = "count", ...) ## S3 method for class 'evinb' terms(x, component = "count", ...)## S3 method for class 'evzinb' logLik(object, ...) ## S3 method for class 'evinb' logLik(object, ...) ## S3 method for class 'evzinb' nobs(object, ...) ## S3 method for class 'evinb' nobs(object, ...) ## S3 method for class 'evzinb' model.frame(formula, ...) ## S3 method for class 'evinb' model.frame(formula, ...) ## S3 method for class 'evzinb' formula(x, component = "count", ...) ## S3 method for class 'evinb' formula(x, component = "count", ...) ## S3 method for class 'evzinb' terms(x, component = "count", ...) ## S3 method for class 'evinb' terms(x, component = "count", ...)
object |
A fitted model. |
... |
Unused. |
formula |
A fitted model (for the |
x |
A fitted model (for the |
component |
Which component ( |
As for the corresponding stats generic. logLik() carries
df = length(object$par.all) and nobs, so AIC() /
BIC() work and match object$AIC / object$BIC.
Fitted values, residuals and simulations from an evzinb / evinb model
## S3 method for class 'evzinb' fitted(object, type = c("harmonic", "explog", "counts", "pareto_alpha"), ...) ## S3 method for class 'evinb' fitted(object, type = c("harmonic", "explog", "counts", "pareto_alpha"), ...) ## S3 method for class 'evzinb' residuals( object, type = c("response", "quantile"), seed = NULL, newdata = NULL, ... ) ## S3 method for class 'evinb' residuals( object, type = c("response", "quantile"), seed = NULL, newdata = NULL, ... ) ## S3 method for class 'evzinb' simulate(object, nsim = 1, seed = NULL, newdata = NULL, ...) ## S3 method for class 'evinb' simulate(object, nsim = 1, seed = NULL, newdata = NULL, ...)## S3 method for class 'evzinb' fitted(object, type = c("harmonic", "explog", "counts", "pareto_alpha"), ...) ## S3 method for class 'evinb' fitted(object, type = c("harmonic", "explog", "counts", "pareto_alpha"), ...) ## S3 method for class 'evzinb' residuals( object, type = c("response", "quantile"), seed = NULL, newdata = NULL, ... ) ## S3 method for class 'evinb' residuals( object, type = c("response", "quantile"), seed = NULL, newdata = NULL, ... ) ## S3 method for class 'evzinb' simulate(object, nsim = 1, seed = NULL, newdata = NULL, ...) ## S3 method for class 'evinb' simulate(object, nsim = 1, seed = NULL, newdata = NULL, ...)
object |
A fitted model. |
type |
For |
... |
For |
seed |
Optional RNG seed for the randomized quantile residuals /
|
newdata |
Optional data to simulate for, or (round10 I.2, for
|
nsim |
Number of simulated response vectors. |
fitted() / residuals() return a numeric vector;
simulate() a data frame with columns sim_1, sim_2,
..., with a "seed" attribute following the
simulate convention (the given seed, or, if
none was given, the RNG state before the draws were made).
Running an extreme value and zero inflated negative binomial model with bootstrapping
evzinb( formula_nb, formula_zi = NULL, formula_evi = NULL, formula_pareto = NULL, data, bootstrap = TRUE, n_bootstraps = 100, multicore = NULL, ncores = NULL, block = NULL, weights = NULL, time = NULL, bootstrap_scheme = NULL, block_length = NULL, family = evinf_family(), boot_seed = NULL, start_seed = NULL, control = evinf_control(), max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.zc.multinomial, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.ZC, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C, verbose = FALSE )evzinb( formula_nb, formula_zi = NULL, formula_evi = NULL, formula_pareto = NULL, data, bootstrap = TRUE, n_bootstraps = 100, multicore = NULL, ncores = NULL, block = NULL, weights = NULL, time = NULL, bootstrap_scheme = NULL, block_length = NULL, family = evinf_family(), boot_seed = NULL, start_seed = NULL, control = evinf_control(), max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.zc.multinomial, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.ZC, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C, verbose = FALSE )
formula_nb |
Formula for the negative binomial (count) component of the model.
May include an |
formula_zi |
Formula for the zero-inflation component of the model. If NULL
taken as the same formula as nb, with any |
formula_evi |
Formula for the extreme-value inflation component of the model.
If NULL taken as the same formula as nb, offset stripped as above. May include
its own |
formula_pareto |
Formula for the pareto (extreme value) component of the
model. If NULL taken as the same formula as nb, offset stripped as above.
|
data |
data to run the model on |
bootstrap |
Should bootstrapping be performed. Needed to obtain standard errors and p-values |
n_bootstraps |
Number of bootstraps to run. For use of bootstrapped p-values, at least 1,000 bootstraps are recommended. For approximate p-values, a lower number can be sufficient |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
block |
Optional case-identifier column for block bootstrapping, given
either as a bare column name ( |
weights |
Optional observation weights (round9 D.2), given as a bare
column name ( |
time |
Optional time index for panel/time-series bootstrap resampling
(round9 F), given as a bare column name ( |
bootstrap_scheme |
One of |
block_length |
Block length for |
family |
An |
boot_seed |
Optional bootstrap seed for reproducibility. When supplied
it is used as-is; when |
start_seed |
Optional seed for the perturbed starts when
|
control |
An |
max.diff.par, max.no.em.steps, max.no.em.steps.warmup, c.lim, prune.c.range, max.upd.par.zc.multinomial, max.upd.par.pl.multinomial, max.upd.par.nb, max.upd.par.pl, no.m.bfgs.steps.multinomial, no.m.bfgs.steps.nb, no.m.bfgs.steps.pl, pdf.pl.type, eta.int, init.Beta.multinom.ZC, init.Beta.multinom.PL, init.Beta.NB, init.Beta.PL, init.Alpha.NB, init.C
|
Deprecated. These EM tuning arguments still work but should be
passed through |
verbose |
Logical: should progress of the full run of the model be tracked? |
An object of class 'evzinb'
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) # Peacekeeping and one-sided violence, with an in-formula log1p() transform # for the Pareto (extreme-value) component (see `?hks`). data(hks) hks_mod <- evzinb( osvAll ~ troopLag + policeLag + militaryobserversLag + epduration + lntpop + brv_AllLag_log + osvAllLagDum + incomp, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE ) summary(hks_mod) glance(hks_mod) predict(hks_mod, type = "harmonic")data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) # Peacekeeping and one-sided violence, with an in-formula log1p() transform # for the Pareto (extreme-value) component (see `?hks`). data(hks) hks_mod <- evzinb( osvAll ~ troopLag + policeLag + militaryobserversLag + epduration + lntpop + brv_AllLag_log + osvAllLagDum + incomp, formula_pareto = ~ log1p(troopLag), data = hks, n_bootstraps = 2, multicore = FALSE ) summary(hks_mod) glance(hks_mod) predict(hks_mod, type = "harmonic")
Stacked tidy / glance for a model comparison
## S3 method for class 'evzinbcomp' tidy(x, ...) ## S3 method for class 'evzinbcomp' glance(x, ...)## S3 method for class 'evzinbcomp' tidy(x, ...) ## S3 method for class 'evzinbcomp' glance(x, ...)
x |
An |
... |
Passed to the per-model |
A tibble with a leading model column.
Bootstrap fits that error are kept as try-error objects rather than
discarded; fits that converge but whose extreme-value tail is so heavy that
summaries from them are effectively unbounded are flagged degenerate
(see alpha_floor in evinf_control). This returns both,
which are useful diagnostics (e.g. a resample with a singular design, or one
whose Pareto shape collapsed).
failed_bootstraps(object)failed_bootstraps(object)
object |
A fitted |
A tibble with columns id, type ("error",
"degenerate" or "not_converged" – round10 G.3: a replicate
that ran without erroring and isn't degenerate, but whose converge
is FALSE, was previously invisible here), message (the
error text, the degeneracy reason, or NA for
"not_converged" – see n_em_steps/c_converged/
c_warmup_capped there instead) and, for "not_converged"
rows, n_em_steps, c_converged, c_warmup_capped
(NA for "error"/"degenerate" rows); zero rows when
every replicate is usable and converged.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) failed_bootstraps(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) failed_bootstraps(model)
A simulated dataset of 1,000 observations with one dependent and three independent variables generated using the EVZINB distribution
genevzinbgenevzinb
## 'genevzinb' A tibble with 1,000 rows and 4 columns:
Dependent variable following the EVZINB distribution
Continuous independent variables following the standard normal distribution
A simulated dataset of 100 observations with one dependent and three independent variables generated using the EVZINB distribution
genevzinb2genevzinb2
## 'genevzinb2' A tibble with 100 rows and 4 columns:
Dependent variable following the EVZINB distribution
Continuous independent variables following the standard normal distribution
EVZINB and EVINB glance functions
## S3 method for class 'evinb' glance(x, ...)## S3 method for class 'evinb' glance(x, ...)
x |
An EVZINB or EVINB object |
... |
Further arguments to be passed to glance() |
A one-row tibble of goodness-of-fit statistics: number of observations
and parameters, family (round9 E.0, e.g. "nbinom/mixture" or "poisson/mixture"),
alpha_NB (NA when the family has none, round9 E.1), C_EV, AIC, BIC,
log-likelihood, whether the EM
algorithm converged (converged) and whether the C_EV profile settled
within max.c.iter (c_converged; NA for a model fitted
before this field existed), the number of EM steps, the number of
observations at or above C_EV, the smallest fitted Pareto shape
(min_alpha_pl; see alpha_pl_floor in
evinf_control), sum_weights (round9 D.2: the sum of
weights =, equal to nobs for an unweighted fit – this,
not nobs, is what aic/bic and the model's degrees of
freedom are computed from), and, for a bootstrapped model, the bootstrap
replicate counts (NA otherwise): n_bootstraps (usable),
n_failed_bootstraps, n_degenerate_bootstraps — these three
partition the number of replicates requested (see evinf_control
for "degenerate") — and n_c_on_boundary, the number of replicates
whose C_EV landed on a candidate-grid endpoint (informational, not counted
as degenerate); oob_fraction_mean/oob_fraction_min/
oob_fraction_max (round10 0.9), the mean/min/max of each usable
replicate's out-of-bag row fraction (NA without bootstraps) –
most informative for the block schemes, where overlapping blocks change
how much of the data a replicate leaves out; n_starts/
n_starts_at_best (round10 G.1, NA for n_starts <= 1),
how many perturbed starts were tried and how many of them reached the
winning log-likelihood (within 1e-4); median_boot_em_steps/
n_boot_c_capped (round10 G.3, NA without bootstraps), the
median EM step count across usable replicates and how many had their
C_EV profile hit max.c.iter without settling in either phase.
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) glance(model)data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) glance(model)
EVZINB and EVINB glance functions
## S3 method for class 'evzinb' glance(x, ...)## S3 method for class 'evzinb' glance(x, ...)
x |
An EVZINB or EVINB object |
... |
Further arguments to be passed to glance() |
A one-row tibble of goodness-of-fit statistics: number of observations
and parameters, family (round9 E.0, e.g. "nbinom/mixture" or "poisson/mixture"),
alpha_NB (NA when the family has none, round9 E.1), C_EV, AIC, BIC,
log-likelihood, whether the EM
algorithm converged (converged) and whether the C_EV profile settled
within max.c.iter (c_converged; NA for a model fitted
before this field existed), the number of EM steps, the number of
observations at or above C_EV, the smallest fitted Pareto shape
(min_alpha_pl; see alpha_pl_floor in
evinf_control), sum_weights (round9 D.2: the sum of
weights =, equal to nobs for an unweighted fit – this,
not nobs, is what aic/bic and the model's degrees of
freedom are computed from), and, for a bootstrapped model, the bootstrap
replicate counts (NA otherwise): n_bootstraps (usable),
n_failed_bootstraps, n_degenerate_bootstraps — these three
partition the number of replicates requested (see evinf_control
for "degenerate") — and n_c_on_boundary, the number of replicates
whose C_EV landed on a candidate-grid endpoint (informational, not counted
as degenerate); oob_fraction_mean/oob_fraction_min/
oob_fraction_max (round10 0.9), the mean/min/max of each usable
replicate's out-of-bag row fraction (NA without bootstraps) –
most informative for the block schemes, where overlapping blocks change
how much of the data a replicate leaves out; n_starts/
n_starts_at_best (round10 G.1, NA for n_starts <= 1),
how many perturbed starts were tried and how many of them reached the
winning log-likelihood (within 1e-4); median_boot_em_steps/
n_boot_c_capped (round10 G.3, NA without bootstraps), the
median EM step count across usable replicates and how many had their
C_EV profile hit max.c.iter without settling in either phase.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) glance(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) glance(model)
zinbboot and nbboot glance functions
## S3 method for class 'zinbboot' glance(x, ...) ## S3 method for class 'nbboot' glance(x, ...)## S3 method for class 'zinbboot' glance(x, ...) ## S3 method for class 'nbboot' glance(x, ...)
x |
An nbboot or zinbboot object |
... |
Further arguments to be passed to glance() |
A one-row tibble of goodness-of-fit statistics, including whether the full-sample model converged and the number of (failed) bootstraps.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) glance(zinb_comp$zinb) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) glance(zinb_comp$nb)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) glance(zinb_comp$zinb) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) glance(zinb_comp$nb)
zipboot and poissonboot glance functions
## S3 method for class 'zipboot' glance(x, ...) ## S3 method for class 'poissonboot' glance(x, ...)## S3 method for class 'zipboot' glance(x, ...) ## S3 method for class 'poissonboot' glance(x, ...)
x |
A poissonboot or zipboot object |
... |
Further arguments to be passed to glance() |
A one-row tibble of goodness-of-fit statistics, including whether the
full-sample model converged and the number of (failed) bootstraps.
alpha is NA: neither baseline has a dispersion parameter
(round9 E.3).
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) glance(zip_comp$zip) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) glance(zip_comp$poisson)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) glance(zip_comp$zip) data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) glance(zip_comp$poisson)
A goodness-of-fit gof tibble for GOF metrics when using modelsummary. The GM tibble can be used to obtain correct table output when making regression tables with modelsummary
gm_evzinbgm_evzinb
## 'gm_evzinb' A tibble with 11 rows and 3 columns:
The modelsummary/broom internal name for the statistic
The table output for the statistic
The number of decimals reported for each statistic by default (can be adapted)
gof_map_evinf(), which builds this object (and lets you
append rows); the regeneration script lives in 'data-raw/gm_evzinb.R'.
modelsummary
Returns the gof_map that modelsummary::modelsummary() needs to
label and format the goodness-of-fit rows produced by glance.evzinb() /
glance.evinb(). The bundled data object gm_evzinb is
simply gof_map_evinf().
gof_map_evinf(extra = NULL)gof_map_evinf(extra = NULL)
extra |
Optional data frame / tibble of additional rows (columns
|
The coefficient table produced by tidy() with
component = "all" has one row per coefficient per component.
modelsummary::modelsummary() cannot render that with its default
arguments (the same limitation it has for nnet::multinom()); pass
shape = term + y.level ~ model together with
gof_map = gof_map_evinf().
A tibble with columns raw (the glance() column name,
unchanged), clean (the title-case row label in the table) and
fmt (default number of decimals). Rows: observations, parameters,
alpha_NB, C_EV, observations at or above C_EV, log-likelihood, AIC, BIC,
the usable / failed / degenerate bootstrap counts, the number of bootstrap
replicates with C_EV on the candidate-grid boundary, and convergence (both
of the EM loop and of the C_EV profile).
gof_map_evinf() gof_map_evinf(extra = data.frame(raw = "n_em_steps", clean = "EM steps", fmt = 0)) data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) if (requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE)) { modelsummary::modelsummary( model, shape = term + y.level ~ model, gof_map = gof_map_evinf() ) }gof_map_evinf() gof_map_evinf(extra = data.frame(raw = "n_em_steps", clean = "EM steps", fmt = 0)) data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) if (requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE)) { modelsummary::modelsummary( model, shape = term + y.level ~ model, gof_map = gof_map_evinf() ) }
A reduced replication data set from Hultman et al. (2013) United Nations Peacekeeping and Civilian Protection in Civil War. Used to reproduce the results of Randahl and Vegelius (2024). To reproduce any other results from Hultman et al. (2013) please download the original replication dataset using the link under source.
hkshks
A tibble with 3746 rows and 13 columns:
The number of observed fatalities from one-sided violence against civilians in the specified conflict-month
The number of UN military troops, in thousands (lagged)
The number of UN police, in thousands (lagged)
The number of UN military observers, in thousands (lagged)
The natural logarithm of the total number of battle-related deaths in the conflict in the previous month, log-transformed as in Hultman et al. (2013)
A dummy variable taking the value 1 if any one-sided violence against civilians took place in the previous conflict month
UCDP/PRIO incompatibility: 1 = territory, 2 = government
The number of months the current conflict-episode has been ongoing
The natural logarithm of the population of the country in which the conflict takes place
'log1p(troopLag)'
'log(epduration)'
'log1p(policeLag)'
'log1p(militaryobserversLag)'
The UN personnel counts ('troopLag', 'policeLag', 'militaryobserversLag') are expressed in **thousands** of personnel. This is a rescaling of the original Hultman, Kathman, and Shannon replication data introduced for Randahl and Vegelius (2024).
The columns ending in '_log' ('troopLag_log', 'policeLag_log', 'militaryobserversLag_log', 'epdur_log', 'brv_AllLag_log') are transforms that were pre-computed for the Pareto component in Randahl and Vegelius (2024). The transform follows the standard rule in the conflict literature: the personnel counts contain structural zeros and use 'log1p()', while an episode duration is always at least one month (at least two in these data) and uses 'log()', the conventional transform for durations. Since evinf 0.9.4 the component formulas accept in-formula transformations, so 'log1p(troopLag)', 'log(epduration)' etc. can be written directly in the formula and these pre-computed columns are no longer needed.
The bundled data carry **no conflict identifier**. The conflict-level cluster
bootstrap reported in Randahl and Vegelius (2024)
therefore cannot be reproduced from the bundled data alone; download the
original replication data (link under 'source') and merge on your own
conflict key to use block in evzinb() / evinb().
https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/6EBCGA
Hultman L, Kathman J, Shannon M (2013).
“United Nations peacekeeping and civilian protection in civil war.”
American Journal of Political Science, 57(4), 875–891.
Randahl D, Vegelius J (2024).
“Inference with Extremes: Accounting for Extreme Values in Count Regression Models.”
International Studies Quarterly, 68(4), sqae137.
doi:10.1093/isq/sqae137.
Likelihood ratio test for individual variables of evzinb
lr_test( object, vars, single = TRUE, bootstrap = FALSE, multicore = NULL, ncores = NULL, exclude_degenerate = TRUE, verbose = FALSE )lr_test( object, vars, single = TRUE, bootstrap = FALSE, multicore = NULL, ncores = NULL, exclude_degenerate = TRUE, verbose = FALSE )
object |
EVZINB or EVINB object to perform likelihood ratio test on |
vars |
Either a list of character vectors with variable names which to be restricted in the LR test or a character vector of variable names. If a list, each character vector of the list will be run separately, allowing for multiple variables to be restricted as once. If a character vector, parameter 'single' can be used to determine whether all variables in the vector should be restricted at once (single = FALSE) or if the variables should be restricted one by one (single = TRUE) |
single |
Logical. Determining whether variables in 'vars' should be restricted individually (single = TRUE) or all at once (single = FALSE) |
bootstrap |
Should LR tests be conducted on each bootstrapped sample or only on the original sample. |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
exclude_degenerate |
Logical. When |
verbose |
Logical. Should the function be verbose? |
The likelihood ratio statistic is
and is compared to a chi-square distribution with degrees of freedom equal to
the number of design-matrix columns removed across all model components (so a
restricted four-level factor contributes three degrees of freedom).
A tibble with one row per performed LR test, or, when bootstrap = TRUE,
a list with the summary tibble (results) and the per-bootstrap statistics
(boot_results).
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) lr_test(model,'x1')data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) lr_test(model,'x1')
Average, over the data, of the numerical effect of each covariate on the chosen predicted quantity, with bootstrap confidence intervals obtained by recomputing the effect on every bootstrap fit.
marginal_effects( object, variables = NULL, type = c("harmonic", "explog", "states", "quantile"), quantile = NULL, at = list(), method = c("derivative", "difference"), eps = 1e-04, delta = 1, conf_level = 0.95, newdata = NULL, n_max = 500, exclude_degenerate = TRUE, multicore = NULL, ncores = NULL, clamp_alpha_pl = FALSE )marginal_effects( object, variables = NULL, type = c("harmonic", "explog", "states", "quantile"), quantile = NULL, at = list(), method = c("derivative", "difference"), eps = 1e-04, delta = 1, conf_level = 0.95, newdata = NULL, n_max = 500, exclude_degenerate = TRUE, multicore = NULL, ncores = NULL, clamp_alpha_pl = FALSE )
object |
A fitted |
variables |
Covariates to compute effects for (default: all raw
covariates). An unrecognised name errors, listing the available
covariates. A variable that only ever enters the model's formulas
wrapped in |
type |
|
quantile |
Quantile for |
at |
Named list of values at which to hold covariates before differencing. |
method |
For numeric covariates, |
eps |
Step size for the central finite difference
( |
delta |
Covariate increase for |
conf_level |
Confidence level for the bootstrap intervals. |
newdata |
Data to average over (default: the estimation data). |
n_max |
For |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate
(default |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
clamp_alpha_pl |
(round11 A4) |
The confidence interval is a percentile interval: conf.low /
conf.high are the empirical quantiles of the effect recomputed on each
bootstrap fit. std.error is the standard deviation of those recomputed
effects, except for type = "quantile" where the bootstrap distribution
is heavy-tailed and std.error is instead a robust scale read off the
interval, (conf.high - conf.low) / (2 * qnorm(1 - (1 - conf_level)/2)).
It is reported for reference and is not used to build the interval.
Harmonic-mean effects and bootstrap intervals. The harmonic-mean
prediction is , which blows up on bootstrap fits
whose Pareto shape is very small. The point estimate for
type = "harmonic" is well behaved, but the percentile interval can be
extremely wide because a handful of bootstrap draws are enormous. For
inference about how a covariate shifts the outcome, prefer
type = "quantile" (effect on a predicted quantile) or
type = "states" (effect on the state probabilities), whose bootstrap
distributions are bounded.
Quantile effects. The mixture quantile is a step function of a
count, so its derivative is either 0 or huge; type = "quantile"
therefore defaults to method = "difference" (the average change in the
predicted quantile from a delta-unit increase in the covariate). The
derivative is still available with method = "derivative" but its
bootstrap distribution can be very heavy-tailed. The quantile effect is
evaluated per observation by bisection on the mixture CDF, which is still
the slowest prediction type, so it is averaged over at most n_max
rows (see that argument).
A tibble with columns variable, contrast, type,
estimate, std.error, conf.low, conf.high.
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 10) marginal_effects(model, variables = "x1")data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 10) marginal_effects(model, variables = "x1")
Out of bag predictive performance of EVZINB and EVINB models
oob_evaluation( object, predict_type = c("harmonic", "explog"), metric = c("rmsle", "rmse", "mse", "mae"), exclude_degenerate = TRUE )oob_evaluation( object, predict_type = c("harmonic", "explog"), metric = c("rmsle", "rmse", "mse", "mae"), exclude_degenerate = TRUE )
object |
A fitted evzinb / evinb model with bootstraps, or an
|
predict_type |
What type of prediction should be made? Harmonic mean, or exp(log(prediction))? |
metric |
What metric should be used for the out of bag evaluation? Default options include rmsle, rmse, mse, and mae. Can also take a user supplied function of the form function(y_pred,y_true) which returns a single value |
exclude_degenerate |
For a single evinf model (or its |
For a single model, a vector of length n_bootstraps. For an
evzinbcomp object, a tibble with one column per compared model
(evinf, nb, zinb, ...) and one row per bootstrap. A
replicate that is NA in any column (a degenerate evinf replicate,
or a compared model's fit that errored) is NA in every
column, so a column-wise na.rm = TRUE summary (e.g.
summarize(oob, across(everything(), median, na.rm = TRUE)))
compares the same set of replicates across models rather than silently
different ones; the row mask itself is available as
attr(out, "excluded") (round8 0.4, review §5).
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) oob_evaluation(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) oob_evaluation(model)
Plot the log-likelihood profile over the candidate values of C_EV
plot_c_profile(object)plot_c_profile(object)
object |
A fitted |
A ggplot object.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) plot_c_profile(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) plot_c_profile(model)
Plots for evzinb / evinb models
## S3 method for class 'evzinb' plot( x, type = c("states", "prediction", "coefficients", "ppc", "ppc_quantiles", "trace"), variable = NULL, quantiles = c(0.5, 0.95), ... ) ## S3 method for class 'evinb' plot( x, type = c("states", "prediction", "coefficients", "ppc", "ppc_quantiles", "trace"), variable = NULL, quantiles = c(0.5, 0.95), ... )## S3 method for class 'evzinb' plot( x, type = c("states", "prediction", "coefficients", "ppc", "ppc_quantiles", "trace"), variable = NULL, quantiles = c(0.5, 0.95), ... ) ## S3 method for class 'evinb' plot( x, type = c("states", "prediction", "coefficients", "ppc", "ppc_quantiles", "trace"), variable = NULL, quantiles = c(0.5, 0.95), ... )
x |
A fitted |
type |
|
variable |
Covariate to vary (required for |
quantiles |
Quantiles to draw for |
... |
Passed to |
A ggplot object.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) plot(model, type = "coefficients") plot(model, type = "prediction", variable = "x1") data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) plot(hks_mod, type = "prediction", variable = "troopLag")data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) plot(model, type = "coefficients") plot(model, type = "prediction", variable = "x1") data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) plot(hks_mod, type = "prediction", variable = "troopLag")
Varies one covariate over a range (or set of values) while holding the others at representative values, and returns the model predictions in long form - the raw material for an effect plot.
predict_grid( object, variable, values = NULL, n = 50, at = list(), fixed = c("mean", "median"), type = c("states", "harmonic", "explog", "quantile", "counts", "pareto_alpha"), quantile = NULL, confint = FALSE, conf_level = 0.9, clamp_alpha_pl = FALSE )predict_grid( object, variable, values = NULL, n = 50, at = list(), fixed = c("mean", "median"), type = c("states", "harmonic", "explog", "quantile", "counts", "pareto_alpha"), quantile = NULL, confint = FALSE, conf_level = 0.9, clamp_alpha_pl = FALSE )
object |
A fitted |
variable |
Name of the covariate to vary. |
values |
Optional explicit values for |
n |
Number of grid points for a numeric |
at |
Named list of values at which to hold specific other covariates. |
fixed |
How to hold the remaining numeric covariates: |
type |
Prediction type: |
quantile |
Quantile for |
confint |
Add bootstrap confidence intervals (not for |
conf_level |
Confidence level. |
clamp_alpha_pl |
(round11 A4) |
A tibble with columns variable, value, type,
estimate, conf.low, conf.high (the last two NA
when confint = FALSE).
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) predict_grid(model, "x1", type = "harmonic") data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) predict_grid(hks_mod, "troopLag", type = "states")data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) predict_grid(model, "x1", type = "harmonic") data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) predict_grid(hks_mod, "troopLag", type = "states")
Predictions from evinb object
## S3 method for class 'evinb' predict( object, newdata = NULL, type = c("harmonic", "explog", "counts", "pareto_alpha", "evinf", "count_state", "states", "all", "quantile", "distribution", "exceedance", "draws"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.95, multicore = NULL, ncores = NULL, return_bootstraps = FALSE, exclude_degenerate = TRUE, support = NULL, max_support = 1e+05, format = c("long", "matrix"), threshold = NULL, n_draws = 1000, parameter_uncertainty = FALSE, seed = NULL, keep = NULL, clamp_alpha_pl = FALSE, ... )## S3 method for class 'evinb' predict( object, newdata = NULL, type = c("harmonic", "explog", "counts", "pareto_alpha", "evinf", "count_state", "states", "all", "quantile", "distribution", "exceedance", "draws"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.95, multicore = NULL, ncores = NULL, return_bootstraps = FALSE, exclude_degenerate = TRUE, support = NULL, max_support = 1e+05, format = c("long", "matrix"), threshold = NULL, n_draws = 1000, parameter_uncertainty = FALSE, seed = NULL, keep = NULL, clamp_alpha_pl = FALSE, ... )
object |
An evinb object for which to produce predicted values |
newdata |
Optional new data (tibble) to produce predicted values from |
type |
Character string, 'harmonic' for the harmonic mean and 'explog' for exponentiated expected log, 'counts' for predicted count of the negative binomial component, 'pareto_alpha' for the predicted pareto alpha value, 'states' for the predicted component states (prior), 'count_state' for predicted probability of the count state, 'evinf' for predicted probability of the pareto state, 'all' for all predicted values, 'quantile' for quantile prediction (scalar or vector, see ‘quantile'), ’distribution' for the full predictive distribution, 'exceedance' for exceedance probabilities (see ‘threshold'), and ’draws' for predictive draws (see 'n_draws'). |
pred |
Type of prediction to be used, defaults to the original prediction from the fitted model, with alternatives being the bootstrapped median or mean. Note that bootstrap mean may yield infinite values, especially when doing quantile prediction |
quantile |
Quantile(s) for which to produce quantile prediction. A single value keeps the existing scalar behavior exactly (a vector, e.g. 'c(.5, .9, .99)' (round10 H.3), returns a tibble with one 'qXX' column per probability ('qXX_lo'/'qXX_hi' too, with 'confint = TRUE')). |
confint |
Should confidence intervals be made for the predictions? Note: only available for vector type predictions and not 'states', 'all' or 'distribution'. |
conf_level |
What confidence level should be used for confidence intervals |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
return_bootstraps |
Should the bootstrapped predictions be returned as well? Useful for further custom analyses of the bootstrapped predictions. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
support |
(round10 H.2) ‘type = ’distribution'' only: shared support (vector of 'y' values) for every row. 'NULL' (default) uses '0:K', 'K' the ceiling of the 0.999 mixture quantile of the heaviest-tailed row, capped at 'max_support'. |
max_support |
(round10 H.2) Upper bound on the default ‘support'’s 'K' (default 1e5); ignored when 'support' is given explicitly. |
format |
(round10 H.2) ‘type = ’distribution'‘ only: '’long'‘ (default; a tibble with '.row', 'y', 'prob') or '’matrix'' (the raw n x K matrix; 'keep' is ignored). |
threshold |
(round10 H.4) ‘type = ’exceedance'' only: a vector of thresholds; result columns are 'p_ge_<threshold>', 'P(Y >= threshold)'. |
n_draws |
(round10 H.5) ‘type = ’draws'' only: number of predictive draws per row. |
parameter_uncertainty |
(round10 H.5) ‘type = ’draws'‘ only: if 'TRUE', each draw uses a randomly chosen usable bootstrap replicate’s parameters instead of the full-sample estimate. |
seed |
(round10 H.5) ‘type = ’draws'‘ only: optional seed; the caller’s '.Random.seed' is left untouched either way (same convention as 'simulate()'). |
keep |
(round10 H.6) Optional character vector of ‘newdata' column names to carry into the result (a join key for panel data); supported for 'type' in '’quantile'‘ (vector), '’distribution'‘ ('format = ’long'‘), '’exceedance'‘ and '’draws''. |
clamp_alpha_pl |
(round11 A4) ‘type = ’explog'‘ (and '’all'‘) only: 'FALSE' (default) uses the fitted 'alpha_pl' as-is and warns below 0.1; 'TRUE' clamps at 0.1; a positive number clamps there instead. Exempt from the global 'alpha_pl_floor' (which stays in force for '’harmonic'‘/'’quantile''). The clamp used is recorded as a '"clamp_alpha_pl"' attribute on the result. |
... |
Not used; any name here is an unknown argument and errors (round11 A5), naming the valid arguments for the requested 'type'. A known argument supplied for a 'type' it does nothing for (e.g. 'threshold' with 'type = "harmonic"') warns instead, since it has a named formal below and never reaches here. |
The likelihood, CDF, quantile prediction, residuals and
simulate() all use the discretised Pareto distribution (the
integer-valued distribution the model is actually fit on). The
'harmonic' and 'explog' point predictions instead use the
harmonic and geometric means of the *continuous* Pareto distribution as a
closed-form approximation to the corresponding discretised-Pareto moments;
this keeps existing point predictions unchanged but means they are not
computed from exactly the same distribution as the rest of the model.
type = 'explog' is , which grows
explosively as the fitted Pareto shape approaches 0.
Unlike 'harmonic' and 'quantile', this path is exempt from
the global evinf_control(alpha_pl_floor = ); clamp_alpha_pl
is its own, separate knob (round11 A4). With clamp_alpha_pl = FALSE
(the default), the fitted is used as-is and a warning
fires whenever any value is below 0.1, since the prediction is
effectively undefined there (, and it
may be Inf); check glance()$min_alpha_pl, and either pass
clamp_alpha_pl = TRUE (clamps at 0.1) or a positive number (clamps
there instead), or use type = 'harmonic'. The clamp actually
applied (FALSE, or the numeric floor) is recorded as a
"clamp_alpha_pl" attribute on the result.
A vector of predicted values for type 'harmonic', 'explog', 'counts', 'pareto_alpha','evinf', 'count_state', and 'quantile' (scalar), or a tibble for type 'states', 'all', 'quantile' (vector), 'distribution' (‘format = ’long'‘), ’exceedance', 'draws', or if confint=T; a matrix for 'distribution' with ‘format = ’matrix''.
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) predict(model) predict(model, type='all', quantile = 0.9) # all available predicted values predict(model, type='quantile', quantile = c(.5, .9, .99)) # several quantiles (H.3) predict(model, type='distribution') # the full predictive distribution (H.2) predict(model, type='exceedance', threshold = c(10, 100)) # exceedance probs (H.4) predict(model, type='draws', n_draws = 100) # predictive draws (H.5)data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) predict(model) predict(model, type='all', quantile = 0.9) # all available predicted values predict(model, type='quantile', quantile = c(.5, .9, .99)) # several quantiles (H.3) predict(model, type='distribution') # the full predictive distribution (H.2) predict(model, type='exceedance', threshold = c(10, 100)) # exceedance probs (H.4) predict(model, type='draws', n_draws = 100) # predictive draws (H.5)
Predictions from evzinb object
## S3 method for class 'evzinb' predict( object, newdata = NULL, type = c("harmonic", "explog", "counts", "pareto_alpha", "zi", "evinf", "count_state", "states", "all", "quantile", "distribution", "exceedance", "draws"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.95, multicore = NULL, ncores = NULL, return_bootstraps = FALSE, exclude_degenerate = TRUE, support = NULL, max_support = 1e+05, format = c("long", "matrix"), threshold = NULL, n_draws = 1000, parameter_uncertainty = FALSE, seed = NULL, keep = NULL, clamp_alpha_pl = FALSE, ... )## S3 method for class 'evzinb' predict( object, newdata = NULL, type = c("harmonic", "explog", "counts", "pareto_alpha", "zi", "evinf", "count_state", "states", "all", "quantile", "distribution", "exceedance", "draws"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.95, multicore = NULL, ncores = NULL, return_bootstraps = FALSE, exclude_degenerate = TRUE, support = NULL, max_support = 1e+05, format = c("long", "matrix"), threshold = NULL, n_draws = 1000, parameter_uncertainty = FALSE, seed = NULL, keep = NULL, clamp_alpha_pl = FALSE, ... )
object |
An evzinb object for which to produce predicted values |
newdata |
Optional new data (tibble) to produce predicted values from |
type |
Character string, 'harmonic' for the harmonic mean and 'explog' for exponentiated expected log, 'counts' for predicted count of the negative binomial component, 'pareto_alpha' for the predicted pareto alpha value, 'states' for the predicted component states (prior), 'count_state' for predicted probability of the count state, 'evinf' for predicted probability of the pareto state,'zi' for the predicted probability of the zero state, 'all' for all predicted values, 'quantile' for quantile prediction (scalar or vector, see ‘quantile'), ’distribution' for the full predictive distribution, 'exceedance' for exceedance probabilities (see ‘threshold'), and ’draws' for predictive draws (see 'n_draws'). |
pred |
Type of prediction to be used, defaults to the original prediction from the fitted model, with alternatives being the bootstrapped median or mean. Note that bootstrap mean may yield infinite values, especially when doing quantile prediction |
quantile |
Quantile(s) for which to produce quantile prediction. A single value keeps the existing scalar behavior exactly (a vector, e.g. 'c(.5, .9, .99)' (round10 H.3), returns a tibble with one 'qXX' column per probability ('qXX_lo'/'qXX_hi' too, with 'confint = TRUE')). |
confint |
Should confidence intervals be made for the predictions? Note: only available for vector type predictions and not 'states', 'all' or 'distribution'. |
conf_level |
What confidence level should be used for confidence intervals |
multicore |
Bootstrap parallelisation shortcut. The default
( |
ncores |
Number of workers when |
return_bootstraps |
Should the bootstrapped predictions be returned as well? Useful for further custom analyses of the bootstrapped predictions. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
support |
(round10 H.2) ‘type = ’distribution'' only: shared support (vector of 'y' values) for every row. 'NULL' (default) uses '0:K', 'K' the ceiling of the 0.999 mixture quantile of the heaviest-tailed row, capped at 'max_support'. |
max_support |
(round10 H.2) Upper bound on the default ‘support'’s 'K' (default 1e5); ignored when 'support' is given explicitly. |
format |
(round10 H.2) ‘type = ’distribution'‘ only: '’long'‘ (default; a tibble with '.row', 'y', 'prob') or '’matrix'' (the raw n x K matrix; 'keep' is ignored). |
threshold |
(round10 H.4) ‘type = ’exceedance'' only: a vector of thresholds; result columns are 'p_ge_<threshold>', 'P(Y >= threshold)'. |
n_draws |
(round10 H.5) ‘type = ’draws'' only: number of predictive draws per row. |
parameter_uncertainty |
(round10 H.5) ‘type = ’draws'‘ only: if 'TRUE', each draw uses a randomly chosen usable bootstrap replicate’s parameters instead of the full-sample estimate. |
seed |
(round10 H.5) ‘type = ’draws'‘ only: optional seed; the caller’s '.Random.seed' is left untouched either way (same convention as 'simulate()'). |
keep |
(round10 H.6) Optional character vector of ‘newdata' column names to carry into the result (a join key for panel data); supported for 'type' in '’quantile'‘ (vector), '’distribution'‘ ('format = ’long'‘), '’exceedance'‘ and '’draws''. |
clamp_alpha_pl |
(round11 A4) ‘type = ’explog'‘ (and '’all'‘) only: 'FALSE' (default) uses the fitted 'alpha_pl' as-is and warns below 0.1; 'TRUE' clamps at 0.1; a positive number clamps there instead. Exempt from the global 'alpha_pl_floor' (which stays in force for '’harmonic'‘/'’quantile''). The clamp used is recorded as a '"clamp_alpha_pl"' attribute on the result. |
... |
Not used; any name here is an unknown argument and errors (round11 A5), naming the valid arguments for the requested 'type'. A known argument supplied for a 'type' it does nothing for (e.g. 'threshold' with 'type = "harmonic"') warns instead, since it has a named formal below and never reaches here. |
The likelihood, CDF, quantile prediction, residuals and
simulate() all use the discretised Pareto distribution (the
integer-valued distribution the model is actually fit on). The
'harmonic' and 'explog' point predictions instead use the
harmonic and geometric means of the *continuous* Pareto distribution as a
closed-form approximation to the corresponding discretised-Pareto moments;
this keeps existing point predictions unchanged but means they are not
computed from exactly the same distribution as the rest of the model.
type = 'explog' is , which grows
explosively as the fitted Pareto shape approaches 0.
Unlike 'harmonic' and 'quantile', this path is exempt from
the global evinf_control(alpha_pl_floor = ); clamp_alpha_pl
is its own, separate knob (round11 A4). With clamp_alpha_pl = FALSE
(the default), the fitted is used as-is and a warning
fires whenever any value is below 0.1, since the prediction is
effectively undefined there (, and it
may be Inf); check glance()$min_alpha_pl, and either pass
clamp_alpha_pl = TRUE (clamps at 0.1) or a positive number (clamps
there instead), or use type = 'harmonic'. The clamp actually
applied (FALSE, or the numeric floor) is recorded as a
"clamp_alpha_pl" attribute on the result.
A vector of predicted values for type 'harmonic', 'explog', 'counts', 'pareto_alpha','zi','evinf', 'count_state', and 'quantile' (scalar), or a tibble for type 'states', 'all', 'quantile' (vector), 'distribution' (‘format = ’long'‘), ’exceedance', 'draws', or if confint=T; a matrix for 'distribution' with ‘format = ’matrix''.
Bootstrap fits (and the per-bootstrap work in add_bootstraps,
lr_test, compare_models,
predict.evzinb and marginal_effects) are
dispatched with furrr on top of a future plan. Set the plan
once for your session and leave multicore at its default:
future::plan(future::multisession, workers = 8) progressr::handlers(global = TRUE) # opt in to progress bars model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 1000)
Any future backend works (multisession, cluster,
callr, a HPC batchtools plan, ...). As a convenience
multicore = TRUE sets a temporary multisession plan for the
single call and restores the previous plan on exit. Large models may need
options(future.globals.maxSize = <bytes>) to raise the default limit
on the data shipped to each worker.
The bootstrap resample indices (object$bootstraps[[i]]$boot_id)
are fully determined by boot_seed and are independent of the
future backend, the number of workers, and chunking — a fixed seed
always draws the same resamples.
The fitted coefficients are reproducible only to within
floating-point noise. BLAS operations are not bit-reproducible across
processes, so a sequential run and a multisession run of the same
seed can differ in the last few digits; and because the EM log-likelihood
can have close local optima, an occasional replicate converges to a
different one under a different backend. For replication material, record
the future plan and utils::sessionInfo() alongside
boot_seed, and produce the canonical set of estimates under a single
fixed plan.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) predict(model) predict(model, type='all', quantile = 0.9) # all available predicted values predict(model, type='quantile', quantile = c(.5, .9, .99)) # several quantiles (H.3) predict(model, type='distribution') # the full predictive distribution (H.2) predict(model, type='exceedance', threshold = c(10, 100)) # exceedance probs (H.4) predict(model, type='draws', n_draws = 100) # predictive draws (H.5)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) predict(model) predict(model, type='all', quantile = 0.9) # all available predicted values predict(model, type='quantile', quantile = c(.5, .9, .99)) # several quantiles (H.3) predict(model, type='distribution') # the full predictive distribution (H.2) predict(model, type='exceedance', threshold = c(10, 100)) # exceedance probs (H.4) predict(model, type='draws', n_draws = 100) # predictive draws (H.5)
Prediction for nbboot
## S3 method for class 'nbboot' predict( object, newdata = NULL, type = c("predicted", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )## S3 method for class 'nbboot' predict( object, newdata = NULL, type = c("predicted", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )
object |
a fitted nbboot object |
newdata |
Data to make predictions on |
type |
What prediction should be computed? |
pred |
Prediction type, 'original', 'bootstrap_median', or 'bootstrap_mean' |
quantile |
Quantile for quantile prediction |
confint |
Should confidence intervals be created? |
conf_level |
Confidence level when predicting with CIs |
... |
Not used |
Predictions from nbboot
Prediction for poissonboot
## S3 method for class 'poissonboot' predict( object, newdata = NULL, type = c("predicted", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )## S3 method for class 'poissonboot' predict( object, newdata = NULL, type = c("predicted", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )
object |
a fitted poissonboot object |
newdata |
Data to make predictions on |
type |
What prediction should be computed? |
pred |
Prediction type, 'original', 'bootstrap_median', or 'bootstrap_mean' |
quantile |
Quantile for quantile prediction |
confint |
Should confidence intervals be created? |
conf_level |
Confidence level when predicting with CIs |
... |
Not used |
Predictions from poissonboot
Prediction for zinbboot
## S3 method for class 'zinbboot' predict( object, newdata = NULL, type = c("predicted", "counts", "zi", "count_state", "states", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )## S3 method for class 'zinbboot' predict( object, newdata = NULL, type = c("predicted", "counts", "zi", "count_state", "states", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )
object |
a fitted zinbboot object |
newdata |
Data to make predictions on |
type |
What prediction should be computed? One of |
pred |
Prediction type, 'original', 'bootstrap_median', or 'bootstrap_mean' |
quantile |
Quantile for quantile prediction |
confint |
Should confidence intervals be created? |
conf_level |
Confidence level when predicting with CIs |
... |
Not used |
Predictions from zinbboot
Prediction for zipboot
## S3 method for class 'zipboot' predict( object, newdata = NULL, type = c("predicted", "counts", "zi", "count_state", "states", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )## S3 method for class 'zipboot' predict( object, newdata = NULL, type = c("predicted", "counts", "zi", "count_state", "states", "all", "quantile"), pred = c("original", "bootstrap_median", "bootstrap_mean"), quantile = NULL, confint = FALSE, conf_level = 0.9, ... )
object |
a fitted zipboot object |
newdata |
Data to make predictions on |
type |
What prediction should be computed? One of |
pred |
Prediction type, 'original', 'bootstrap_median', or 'bootstrap_mean' |
quantile |
Quantile for quantile prediction |
confint |
Should confidence intervals be created? |
conf_level |
Confidence level when predicting with CIs |
... |
Not used |
Predictions from zipboot
EVINB print function
## S3 method for class 'evinb' print(x, ...)## S3 method for class 'evinb' print(x, ...)
x |
A fitted evinb model |
... |
Not used |
x, invisibly.
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(model)data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(model)
EVZINB print function
## S3 method for class 'evzinb' print(x, ...)## S3 method for class 'evzinb' print(x, ...)
x |
A fitted evzinb model |
... |
Not used |
x, invisibly.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(model)
Shows the compared models, a compact glance() table (one row
per slot, key fit statistics only), and the paired bootstrap comparison
from compare_fit() (round10 I.4, audit §5.9) – previously
only the compared-model names and bootstrap count.
## S3 method for class 'evzinbcomp' print(x, metrics = c("aic", "bic", "rmse", "rmsle"), ...)## S3 method for class 'evzinbcomp' print(x, metrics = c("aic", "bic", "rmse", "rmsle"), ...)
x |
An |
metrics |
Metrics for the |
... |
Not used |
x, invisibly.
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(compare_models(model))data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) print(compare_models(model))
Print methods for evzinb / evinb summaries
## S3 method for class 'summary.evzinb' print( x, digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), ... ) ## S3 method for class 'summary.evinb' print( x, digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), ... )## S3 method for class 'summary.evzinb' print( x, digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), ... ) ## S3 method for class 'summary.evinb' print( x, digits = max(3L, getOption("digits") - 3L), signif.stars = getOption("show.signif.stars"), ... )
x |
A |
digits |
Number of significant digits for the coefficient tables. |
signif.stars |
Logical; show significance stars. |
... |
Not used. |
x, invisibly.
Superseded: kept for backwards compatibility, with no plans for
removal, but new code should use simulate(),
which returns a tidy data frame (and a reproducible "seed" attribute)
instead of a bare vector/list.
revinb_fit(object, newdata = NULL, n_draws = 1)revinb_fit(object, newdata = NULL, n_draws = 1)
object |
A fitted EVINB object |
newdata |
Optional newdata |
n_draws |
Number of random draws to make |
A vector of randomly drawn values from the fitted evinb if n_draws == 1, or a list of length n_draws with random drawn values if n_draws > 1
data(genevzinb2) model <- evinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) revinb_fit(model)data(genevzinb2) model <- evinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) revinb_fit(model)
Superseded: kept for backwards compatibility, with no plans for
removal, but new code should use simulate(),
which returns a tidy data frame (and a reproducible "seed" attribute)
instead of a bare vector/list.
revzinb_fit(object, newdata = NULL, n_draws = 1)revzinb_fit(object, newdata = NULL, n_draws = 1)
object |
A fitted EVZINB object |
newdata |
Optional newdata |
n_draws |
Number of random draws to make |
A vector of randomly drawn values from the fitted evzinb if n_draws == 1, or a list of length n_draws with random drawn values if n_draws > 1
data(genevzinb2) model <- evzinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) revzinb_fit(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3, data=genevzinb2, n_bootstraps = 5) revzinb_fit(model)
Cross-tabulate prior vs posterior state classification
state_table(object)state_table(object)
object |
A fitted |
An object of class evinf_state_table: the map_prior x
map_posterior contingency table with row percentages.
data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) state_table(model)data(genevzinb2) model <- evzinb(y ~ x1 + x2 + x3, data = genevzinb2, n_bootstraps = 5) state_table(model)
EVINB summary function
## S3 method for class 'evinb' summary( object, coef = c("original", "bootstrapped_mean", "bootstrapped_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "both", "none"), bootstrapped_props = c("none", "mean", "median"), approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )## S3 method for class 'evinb' summary( object, coef = c("original", "bootstrapped_mean", "bootstrapped_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "both", "none"), bootstrapped_props = c("none", "mean", "median"), approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )
object |
an EVINB object |
coef |
Type of coefficients. Original are the coefficient estimates from the non-bootstrapped version of the model. 'bootstrapped_mean' are the mean coefficients across bootstraps, and 'bootstrapped_median' are the median coefficients across bootstraps |
standard_error |
Should standard errors be computed? |
p_value |
What type of p_values should be computed? 'bootstrapped' are bootstrapped p_values through confidence interval inversion. 'approx' are p-values based on the t-value produced by dividing the coefficient with the standard error. 'both' returns both. |
bootstrapped_props |
Type of bootstrapped proportions of component proportions to be returned |
approx_t_value |
Should approximate t-values be returned |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
... |
Additional arguments passed to the summary function |
When object was fitted with bootstrap = FALSE, the
bootstrap-based quantities (standard errors, p-values, approximate t-values
and bootstrapped proportions) are unavailable; the summary then reports the
point estimates only, n_failed_bootstraps is NA, and a message
is emitted. Requesting bootstrapped coefficients for such a model is an error.
A bootstrapped p-value is never reported below 1 / B, where
B is the number of usable bootstrap replicates: with no draw
crossing the estimate, the true p-value could be anywhere in
[0, 1/B), so it is floored at 1/B rather than reported as
exactly 0. print.summary.evzinb() / print.summary.evinb()
show this as "< 1/B" (e.g. "<0.01" for 100 usable
bootstraps) rather than the default, misleadingly precise "<2e-16".
An EVINB summary object
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) summary(model)data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) summary(model)
EVZINB summary function
## S3 method for class 'evzinb' summary( object, coef = c("original", "bootstrapped_mean", "bootstrapped_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "both", "none"), bootstrapped_props = c("none", "mean", "median"), approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )## S3 method for class 'evzinb' summary( object, coef = c("original", "bootstrapped_mean", "bootstrapped_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "both", "none"), bootstrapped_props = c("none", "mean", "median"), approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )
object |
an EVZINB object |
coef |
Type of coefficients. Original are the coefficient estimates from the non-bootstrapped version of the model. 'bootstrapped_mean' are the mean coefficients across bootstraps, and 'bootstrapped_median' are the median coefficients across bootstraps |
standard_error |
Should standard errors be computed? |
p_value |
What type of p_values should be computed? 'bootstrapped' are bootstrapped p_values through confidence interval inversion. 'approx' are p-values based on the t-value produced by dividing the coefficient with the standard error. 'both' returns both. |
bootstrapped_props |
Type of bootstrapped proportions of component proportions to be returned |
approx_t_value |
Should approximate t-values be returned |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
... |
Additional arguments passed to the summary function |
When object was fitted with bootstrap = FALSE, the
bootstrap-based quantities (standard errors, p-values, approximate t-values
and bootstrapped proportions) are unavailable; the summary then reports the
point estimates only, n_failed_bootstraps is NA, and a message
is emitted. Requesting bootstrapped coefficients for such a model is an error.
A bootstrapped p-value is never reported below 1 / B, where
B is the number of usable bootstrap replicates: with no draw
crossing the estimate, the true p-value could be anywhere in
[0, 1/B), so it is floored at 1/B rather than reported as
exactly 0. print.summary.evzinb() / print.summary.evinb()
show this as "< 1/B" (e.g. "<0.01" for 100 usable
bootstraps) rather than the default, misleadingly precise "<2e-16".
An EVZINB summary object
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) summary(model)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) summary(model)
EVINB tidy function
## S3 method for class 'evinb' tidy( x, component = c("all", "count", "evi", "pareto"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )## S3 method for class 'evinb' tidy( x, component = c("all", "count", "evi", "pareto"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )
x |
An evinb object |
component |
Which component should be shown? One of |
coef_type |
Type of coefficients. Original are the coefficient estimates from the non-bootstrapped version of the model. 'bootstrap_mean' are the mean coefficients across bootstraps, and 'bootstrap_median' are the median coefficients across bootstraps |
standard_error |
Should bootstrapped standard errors be computed? |
p_value |
What type of p_values should be computed? 'bootstrapped' are bootstrapped p_values through confidence interval inversion. 'approx' are p-values based on the t-value produced by dividing the coefficient with the standard error. |
confint |
What type of confidence interval should be computed: |
conf_level |
Confidence level for the confidence interval |
approx_t_value |
Should approximate t-values be returned |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
... |
Other arguments passed to the tidy function |
When x was fitted with bootstrap = FALSE, standard
errors, p-values, t-values and confidence intervals are unavailable; the
estimate column is returned on its own and a message is emitted. Requesting
bootstrapped coefficients for such a model is an error.
A bootstrapped p.value is never reported below 1 / B, where
B is the number of usable bootstrap replicates: with no draw
crossing the estimate, the true p-value could be anywhere in
[0, 1/B), so it is floored at 1/B rather than reported as
exactly 0.
With component = "all" the output has one row per coefficient
per component (a y.level column names the component). Like
nnet::multinom(), this multi-component shape cannot be rendered by
modelsummary::modelsummary() with its default arguments; pass
shape = term + y.level ~ model so each component becomes its own
block of rows.
A tibble with one row per coefficient
data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) tidy(model)data(genevzinb2) model <- evinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) tidy(model)
EVZINB tidy function
## S3 method for class 'evzinb' tidy( x, component = c("all", "count", "zero", "evi", "pareto"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )## S3 method for class 'evzinb' tidy( x, component = c("all", "count", "zero", "evi", "pareto"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, exclude_degenerate = TRUE, ... )
x |
An evzinb object |
component |
Which component should be shown? One of |
coef_type |
Type of coefficients. Original are the coefficient estimates from the non-bootstrapped version of the model. 'bootstrap_mean' are the mean coefficients across bootstraps, and 'bootstrap_median' are the median coefficients across bootstraps |
standard_error |
Should bootstrapped standard errors be computed? |
p_value |
What type of p_values should be computed? 'bootstrapped' are bootstrapped p_values through confidence interval inversion. 'approx' are p-values based on the t-value produced by dividing the coefficient with the standard error. |
confint |
What type of confidence interval should be computed: |
conf_level |
Confidence level for the confidence interval |
approx_t_value |
Should approximate t-values be returned |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate (default TRUE); see the alpha_floor argument of evinf_control(). |
... |
Other arguments passed to the tidy function |
When x was fitted with bootstrap = FALSE, standard
errors, p-values, t-values and confidence intervals are unavailable; the
estimate column is returned on its own and a message is emitted. Requesting
bootstrapped coefficients for such a model is an error.
A bootstrapped p.value is never reported below 1 / B, where
B is the number of usable bootstrap replicates: with no draw
crossing the estimate, the true p-value could be anywhere in
[0, 1/B), so it is floored at 1/B rather than reported as
exactly 0.
With component = "all" the output has one row per coefficient
per component (a y.level column names the component). Like
nnet::multinom(), this multi-component shape cannot be rendered by
modelsummary::modelsummary() with its default arguments; pass
shape = term + y.level ~ model so each component becomes its own
block of rows.
A tibble with one row per coefficient
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) tidy(model) # multi-component table with modelsummary if (requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE)) { modelsummary::modelsummary( model, shape = term + y.level ~ model, gof_map = gof_map_evinf() ) } data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) tidy(hks_mod)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) tidy(model) # multi-component table with modelsummary if (requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE)) { modelsummary::modelsummary( model, shape = term + y.level ~ model, gof_map = gof_map_evinf() ) } data(hks) hks_mod <- evzinb(osvAll ~ troopLag + lntpop + brv_AllLag_log, data = hks, n_bootstraps = 2, multicore = FALSE) tidy(hks_mod)
Tidy function for nbboot
## S3 method for class 'nbboot' tidy( x, coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, include_ylev = FALSE, ... )## S3 method for class 'nbboot' tidy( x, coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, include_ylev = FALSE, ... )
x |
A fitted bootstrapped zero-inflated model |
coef_type |
What type of coefficient should be reported, original, bootstrapped mean, or bootstrapped median |
standard_error |
Should bootstrapped standard errors be reported? |
p_value |
What type of p-value should be reported? Bootstrapped p_values, approximate p-values, or none? |
confint |
What type of confidence intervals should be reported? Bootstrapped p_values, approximate p-values, or none? |
conf_level |
Confidence level for confidence intervals |
approx_t_value |
Should approximate t_values be reported |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
include_ylev |
Logical. Should y.lev be included in the tidy output? Makes for nicer tables when using modelsummary |
... |
Other arguments to be passed to tidy |
A tidy function for a bootstrapped nb model
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) tidy(zinb_comp$nb)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) tidy(zinb_comp$nb)
Tidy function for poissonboot
## S3 method for class 'poissonboot' tidy( x, coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, include_ylev = FALSE, ... )## S3 method for class 'poissonboot' tidy( x, coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, include_ylev = FALSE, ... )
x |
A fitted bootstrapped Poisson model |
coef_type |
What type of coefficient should be reported, original, bootstrapped mean, or bootstrapped median |
standard_error |
Should bootstrapped standard errors be reported? |
p_value |
What type of p-value should be reported? Bootstrapped p_values, approximate p-values, or none? |
confint |
What type of confidence intervals should be reported? Bootstrapped p_values, approximate p-values, or none? |
conf_level |
Confidence level for confidence intervals |
approx_t_value |
Should approximate t_values be reported |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
include_ylev |
Logical. Should y.lev be included in the tidy output? Makes for nicer tables when using modelsummary |
... |
Other arguments to be passed to tidy |
A tidy function for a bootstrapped Poisson model
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) tidy(zip_comp$poisson)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) tidy(zip_comp$poisson)
Tidy function for zinbboot
## S3 method for class 'zinbboot' tidy( x, component = c("all", "count", "zero"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, ... )## S3 method for class 'zinbboot' tidy( x, component = c("all", "count", "zero"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, ... )
x |
A fitted bootstrapped zero-inflated model |
component |
Which component should be shown? |
coef_type |
What type of coefficient should be reported, original, bootstrapped mean, or bootstrapped median |
standard_error |
Should bootstrapped standard errors be reported? |
p_value |
What type of p-value should be reported? Bootstrapped p_values, approximate p-values, or none? |
confint |
What type of confidence intervals should be reported? Bootstrapped p_values, approximate p-values, or none? |
conf_level |
Confidence level for confidence intervals |
approx_t_value |
Should approximate t_values be reported |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
... |
Other arguments to be passed to tidy |
A tidy function for a bootstrapped zinb model
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) tidy(zinb_comp$zinb)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, n_bootstraps = 5) zinb_comp <- compare_models(model) tidy(zinb_comp$zinb)
Tidy function for zipboot
## S3 method for class 'zipboot' tidy( x, component = c("all", "count", "zero"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, ... )## S3 method for class 'zipboot' tidy( x, component = c("all", "count", "zero"), coef_type = c("original", "bootstrap_mean", "bootstrap_median"), standard_error = TRUE, p_value = c("bootstrapped", "approx", "none"), confint = c("none", "bootstrapped", "approx"), conf_level = 0.95, approx_t_value = TRUE, symmetric_bootstrap_p = TRUE, ... )
x |
A fitted bootstrapped zero-inflated Poisson model |
component |
Which component should be shown? |
coef_type |
What type of coefficient should be reported, original, bootstrapped mean, or bootstrapped median |
standard_error |
Should bootstrapped standard errors be reported? |
p_value |
What type of p-value should be reported? Bootstrapped p_values, approximate p-values, or none? |
confint |
What type of confidence intervals should be reported? Bootstrapped p_values, approximate p-values, or none? |
conf_level |
Confidence level for confidence intervals |
approx_t_value |
Should approximate t_values be reported |
symmetric_bootstrap_p |
Should bootstrap p-values be computed as symmetric (leaving alpha/2 percent in each tail)? FALSE gives non-symmetric, but narrower, intervals. TRUE corresponds most closely to conventional p-values. |
... |
Other arguments to be passed to tidy |
A tidy function for a bootstrapped ZIP model
data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) tidy(zip_comp$zip)data(genevzinb2) model <- evzinb(y~x1+x2+x3,data=genevzinb2, family = "poisson", n_bootstraps = 5) zip_comp <- compare_models(model) tidy(zip_comp$zip)
Update and re-fit an evzinb / evinb model
## S3 method for class 'evzinb' update( object, formula_nb., formula_zi., formula_evi., formula_pareto., ..., evaluate = TRUE ) ## S3 method for class 'evinb' update( object, formula_nb., formula_evi., formula_pareto., ..., evaluate = TRUE )## S3 method for class 'evzinb' update( object, formula_nb., formula_zi., formula_evi., formula_pareto., ..., evaluate = TRUE ) ## S3 method for class 'evinb' update( object, formula_nb., formula_evi., formula_pareto., ..., evaluate = TRUE )
object |
A fitted model (must carry |
formula_nb., formula_zi., formula_evi., formula_pareto.
|
Optional
|
... |
|
evaluate |
If |
If object's candidate range for was itself
data-driven (control$c.lim was NULL at the original fit;
see object$c_lim_default), and data is one of the arguments
being changed, control$c.lim and control$init.C are reset to
NULL so they are re-resolved from the new data (with the usual
message) instead of silently reusing the range chosen for the old data. A
c.lim the user pinned explicitly is always kept as-is.
The updated fit, or the call.
Bootstrap covariance matrix of an evzinb / evinb model
## S3 method for class 'evzinb' vcov(object, exclude_degenerate = TRUE, ...) ## S3 method for class 'evinb' vcov(object, exclude_degenerate = TRUE, ...)## S3 method for class 'evzinb' vcov(object, exclude_degenerate = TRUE, ...) ## S3 method for class 'evinb' vcov(object, exclude_degenerate = TRUE, ...)
object |
A fitted model with bootstraps. |
exclude_degenerate |
Drop bootstrap replicates flagged degenerate
(default |
... |
Unused. |
c_ev is estimated on the grid of unique observed values of
the response, not by a smooth optimiser, so it sits on a discrete scale.
Its bootstrap variance (and a percentile confint() on it) is
meaningful, but delta-method standard errors that perturb c_ev
continuously – as marginaleffects::avg_slopes() etc. do – should
be read with that in mind. For uncertainty on covariate effects,
prefer the bootstrap route, marginal_effects().
The covariance matrix of the bootstrap coefficient draws, with the
same names and order as coef(object) (includes alpha_nb and
c_ev).