This is the first CRAN release since 0.8.10 (2024). 0.9.x, 0.10.0 and 0.11.0
were development versions that were never released to CRAN; their NEWS
sections are kept below, unmerged and unrewritten, as the development record.
Everything in this section is therefore, as far as CRAN is concerned, a
change since 0.8.10 -- see cran-comments.md for that framing. This round
(round 11) closes the five release blockers found in a post-round-10 CRAN-
readiness review and prepares the package for submission.
predict() (and fitted(), predict_grid(), marginal_effects(), and
the marginaleffects get_predict() method) now validate ...: an
argument name that is not a known predict() argument at all errors,
naming the valid arguments for the requested type; a known argument
that does nothing for the requested type (e.g. threshold with
type = "harmonic") warns instead of being silently ignored. Previously
a typo, or an argument passed for the wrong type, was accepted silently
and produced a default-argument result (round11 A5).predict(confint = TRUE) (and pred = "bootstrap_median"/
"bootstrap_mean") now errors clearly when every bootstrap replicate is
unusable, instead of crashing with a confusing dplyr "column id not
found" error (found while testing A4 on hks).marginaleffects (>= 0.22.0) in Suggests and guarded the one
affected test: avg_predictions()/predictions() with a model-specific
type ("quantile", "states", "exceedance") errors on
marginaleffects older than 0.22.0, since older versions validate type
against an internal dictionary that does not exempt an unregistered model
class the way 0.22.0 on does (round11 A3).VignetteBuilder: quarto, the quarto entry in Suggests, and
Config/Needs/website: quarto until the JSS vignette actually ships --
R CMD build failed outright without quarto installed, and
R CMD check --as-cran reported a VignetteBuilder-without-a-prebuilt-
index NOTE (round11 A2).\dontrun{} to \donttest{}
(?evzinb, ?evinb, ?plot.evzinb, ?predict_grid, ?tidy.evzinb,
?compare_models), trimming n_bootstraps so --run-donttest stays
fast (round11 B1).summary() no longer errors ("Repaired names have length 3 instead of
length 1") when exactly one bootstrap replicate is usable -- found via
--run-donttest on ?evzinb's own example after the previous fix
unhid it (round11 B1 follow-up).predict()/fitted() gain clamp_alpha_pl for type = "explog"
(also reachable through type = "all", predict_grid() and
marginal_effects()): FALSE (the default) uses the fitted alpha_pl
as-is, exempt from the global alpha_pl_floor (which stays in force for
"harmonic"/"quantile"), and warns when any value is below 0.1, since
C * exp(1 / alpha_pl) is effectively undefined there and may be Inf;
TRUE clamps at 0.1; a positive number clamps there instead. The clamp
actually used is recorded as a "clamp_alpha_pl" attribute on the
result. Previously the explog path silently ran through the (much
tighter) alpha_pl_floor, producing an astronomically large but finite
value instead of Inf (round11 A4).Implements audit §5.6 (offsets and weights), §5.5 (model families) and §5.7
(panel/time-series resampling): work packages D, E and F of
dev/plan_0.12_families_offsets_panel.md.
object$fitted$y.hat.pl_exp.E.logy, y.hat.pl_E.inv.y,
y.hat.pl_median and y.hat.pl_mean fields are no longer stored. They
duplicated predict(type = "explog"/"harmonic"/"quantile") but were not
hurdle-aware (on a hurdle fit, correlation with the corresponding
predict() output was as low as 0.9988). fitted() already routes
through predict(); use that instead. The other object$fitted$* fields
(mu.nb, alpha.pl, the pl_*/prob_*/posterior_* summaries) are
unaffected (round10 0.6).Round-9 follow-ups (dev/review_round9.md):
lr_test() refits now carry the full model's weights, family and
every offset() term. Previously a weighted fit's restricted
log-likelihood came from an unweighted refit (silently wrong, occasionally
a negative LR statistic), a Poisson or hurdle fit's restricted refit
errored or compared against the wrong nested model, and an offset model's
restricted formula silently dropped its offset (round10 0.1).compare_models()'s NB/ZINB/Poisson/ZIP competitor baselines -- including
winsorised/razorised variants and every bootstrap refit -- are now fitted
with the same weights as the evzinb/evinb model, so compare_fit()'s
AIC/BIC comparisons are no longer between a weighted and an unweighted
likelihood (round10 0.2).evzinb()/evinb()/add_bootstraps() now validate time against
block once, up front, before any fitting. Previously a duplicated or
out-of-order time under bootstrap_scheme = "moving_block"/
"stationary" let the full-sample fit complete and only surfaced as every
bootstrap replicate failing; the error now fires immediately and, when
time was given without block, names that as the likely cause
(round10 0.4).b63b302 for genevzinb2 and
hks, both evzinb() and evinb()); a residual of at most a couple of
ULPs can remain on other inputs, down from up to 8.6e-8. See
has_weights in src/evinf.cpp and test-weights.R's round10 0.5 tests
(round10 0.5).predict(type = "explog") warns when any fitted Pareto alpha_pl used in
the prediction is below 0.1: the geometric-mean prediction
C * exp(1 / alpha_pl) is effectively undefined well before that (and
before evinf_control(alpha_pl_floor = ), default 0.01, clamps it to a
finite but still absurd value). No separate clamp is applied; see
?predict.evzinb (round10 0.6).nbinom) fit no longer emits a stream of
NA/Inf replaced by maximum positive value warnings from the NB
dispersion line search. The search's eta interval is now shrunk (with a
small margin) so it never probes a negative alpha_nb, instead of relying
on stats::optimise() to silently recover from evaluating the
log-likelihood there; default-family fits are unaffected (round10 0.7).glance() reports oob_fraction_mean/oob_fraction_min/
oob_fraction_max (NA without bootstraps), a fit-level summary across
usable bootstrap replicates of a quantity that previously lived only on
each replicate; most informative for the block schemes, where overlapping
blocks change how much of the data a replicate leaves out (round10 0.9).time given without block, with repeated values, now warns (rather
than passing silently) for bootstrap_schemes that don't use time
("iid"/"cluster") -- "moving_block"/"stationary" already error on
this (round10 0.4); time is otherwise unused for the other schemes, so
this is the same likely-missing-block = mistake surfacing outside the
case that errors (round10 0.9).pscl sign
convention (hurdle()'s zero component models P(Y > 0), evinf's models
the zero state directly) in ?evinf_family (round10 0.9).R-CMD-check.yaml's dependency-install step now retries once on
failure and pins an explicit cache key, to absorb the RSPM binary
download that has intermittently failed on its first attempt (seen on
macOS). Added ^\.claude$ to .Rbuildignore so a local .claude
directory cannot produce an R CMD check NOTE or leak into a submitted
tarball. test-coverage.yaml already references
secrets.CODECOV_TOKEN; that repository secret still needs to be created
in GitHub settings for Codecov to report a figure (round10 0.8).evzinb() / evinb() project total bootstrap runtime from the first
min(4, n_bootstraps) replicates actually completing (round10 J.1, audit
§5.11), replacing the old "full-sample fit time x n_bootstraps" proxy,
which ignored parallelism entirely and could differ substantially from a
bootstrap replicate's own fit time. Shown (via message()) when
verbose = TRUE or the projection exceeds 60 seconds. The probe
replicates duplicate (rather than change) part of the real dispatch's own
work, so bootstrap output for a given boot_seed is unaffected.
evinf_control(chunk_size = ) (round10 J.2, audit §5.11): the number of
replicates/starts evinf_pmap() hands a worker at a time. NULL (the new
default, replacing a hardcoded 1L) picks
max(1, ceiling(B / (4 * future::nbrOfWorkers()))) at dispatch time --
four chunks per worker. Chunk size only changes how work is grouped for
dispatch, never each replicate's RNG stream (verified empirically), so
bootstrap/start output is identical across chunk sizes for the same seed.
evinf_bench_data(n, seed) (round10 J.3, audit §5.11, decision D8): a
synthetic-data generator with known true EVZINB parameters, for
benchmarking and scale testing at whatever n is needed without bundling
an additional .rda per size. inst/bench/bench_evinf.R now uses it
(replacing the bundled genevzinb2/hks data it previously benchmarked
against) and compares each timing to a stored reference
(inst/bench/reference_timings.rds), flagging (in the CI log only --
never failing the build) a case more than 2x slower than its reference.
marginaleffects support (round10 I.1, audit §5.9) gains type = "states" (a long group/estimate frame, one group per prior state --
evinb has no zero state), "quantile" (a single quantile = passed
through ..., defaulting to the median) and "exceedance" (a single or
several threshold = passed through ..., fanning out into groups like
"states" when more than one is given). type = "quantile" uses the
continuous (linearly-interpolated) mixture quantile for delta-method
standard errors, the same fix marginal_effects()'s own
predict_from_boot() already applied (audit N1) -- the integer quantile's
derivative is zero almost everywhere, which otherwise starves the delta
method into NA. marginaleffects does not know evinf's extra ...
arguments and will warn "not known to be supported" for quantile =/
threshold =; the prediction is still computed correctly.
augment.evzinb() / augment.evinb() (round10 I.2, audit §5.9): a
generics::augment()/broom-style method appending .fitted
(predict(type = "harmonic")), .prob_zero/.prob_count/.prob_evi
(prior state probabilities; evinb has no zero state), .state (the MAP
prior state), and, when the response is present, .resid (a seeded
randomized quantile residual), .post_zero/.post_count/.post_evi
(posterior state probabilities) and .post_state (the MAP posterior
state). residuals(type = ) gains a newdata = argument to support
this (previously always the estimation data; unchanged when omitted).
print.evzinbcomp() (round10 I.4, audit §5.9) now also shows a compact
glance() table (one row per compared model: nobs, npar, logLik,
aic, bic) and the compare_fit() paired-bootstrap comparison --
previously only the compared-model names and bootstrap count.
Internal evinf_pmf() / evinf_cdf() (round10 H.1, audit §5.8) are now
the package's one definition of the predictive distribution, family-aware
across all four count x zero combinations. mixture_p()'s CDF formula and
a new shared per-state-density building block (evinf_state_densities())
back both; classify_states(newdata = ) and residuals(type = "quantile") are refactored onto them, and as a result now also clamp a
collapsed alpha_pl against alpha_pl_floor (with the same warning as
every predict() type already gives) -- they did not before. This sets up
predict(type = "distribution"/"quantile"/"exceedance"/"draws") (below).
evinf_dist_params()/evinf_pmf()/evinf_cdf() also gain an explicit
is_zinb = override (default inherits(object, "evzinb"), correct for a
full-sample fit): a bootstrap replicate does not reliably carry that
class, so the default silently misrouted an evzinb replicate to the
evinb parameter extraction. Every caller in this package that runs one
of these on object$bootstraps now passes is_zinb explicitly, the same
pattern the existing bootstrap-prediction code already used.
predict(type = "distribution") (round10 H.2): the full predictive
distribution as a long tibble (.row, y, prob; format = "matrix"
for the raw n x K matrix). support defaults to 0:K, K the ceiling of
the 0.999 mixture quantile of the heaviest-tailed row, capped at
max_support (default 1e5). No confint support (the distribution is
already the quantity a CI would summarize).
predict(type = "quantile", quantile = <vector>) (round10 H.3): several
mixture quantiles at once (e.g. quantile = c(.5, .9, .99)), returned as
a tibble with one qXX column per probability (qXX_lo/qXX_hi too,
with confint = TRUE), sharing one bisection over all probabilities at
once rather than one call per probability. A length-1 quantile is
unaffected -- still the original scalar return.
predict(type = "exceedance", threshold = <vector>) (round10 H.4):
exceedance probabilities P(Y >= threshold), as p_ge_<threshold>
columns (via 1 - F(threshold - 1); threshold = 0 is exactly 1 for
every count-valued model here, computed directly rather than at the
undefined F(-1)). confint = TRUE supported via the same bootstrap
machinery as the existing vector-output types.
predict(type = "draws", n_draws =, parameter_uncertainty =) (round10
H.5): predictive draws as a long tibble (.row, .draw, y), sharing
the sampler simulate() already uses (revzinb_fit()/revinb_fit())
rather than a second one. parameter_uncertainty = TRUE draws each
observation's parameters from a randomly chosen usable bootstrap
replicate instead of the full-sample estimate. Seeded like simulate();
leaves the caller's .Random.seed untouched either way.
keep = c(...) (round10 H.6): every H.2-H.5 type above accepts a
character vector of newdata column names to carry into the result, for
a join key on panel data instead of relying on row order.
evzinb() / evinb() gain evinf_control(n_starts =, start_jitter =)
(round10 G.1, audit §5.10): with n_starts > 1, the default start plus
n_starts - 1 perturbed starts (coefficients jittered by
N(0, start_jitter^2), C_EV drawn uniformly from the candidate grid)
each run the EM in parallel (via the same multicore =/ncores = the
bootstrap dispatch uses), and the highest-log-likelihood replicate is
kept -- the EM can land in different local optima depending on
floating-point accumulation order, and multiple starts are the principled
response. Every replicate's outcome is recorded in object$starts
(NULL for the default n_starts = 1, to keep single-start objects the
same size as before); print()/glance() report how many starts reached
the winning log-likelihood. Bootstrap replicates warm-start from the
winning solution, same as a single-start fit. A new start_seed =
argument on evzinb()/evinb() seeds the perturbed starts and is
recorded as object$start_seed.
plot(type = "trace") (round10 G.2): the log-likelihood and C_EV
traces, in two stacked panels with the warm-up/convergence boundary
marked. With n_starts > 1, every perturbed start's trace is overlaid in
grey behind the winning start's, in black.
Per-bootstrap convergence detail (round10 G.3, audit §5.10): every
bootstrap replicate now also carries n_em_steps (converge/
c_converged/c_warmup_capped were already there). failed_bootstraps()
gains a "not_converged" category for a replicate that ran without
erroring and isn't degenerate, but didn't converge -- previously
invisible there. glance() gains median_boot_em_steps and
n_boot_c_capped (replicates whose C_EV profile hit max.c.iter
without settling, in either phase).
check_evinf() (round10 G.4, audit §5.10): a single diagnostic report
(a classed tibble with a print() method) gathering convergence,
C_EV-on-a-boundary, the smallest fitted Pareto shape against
alpha_pl_floor, failed/degenerate bootstraps, start agreement
(n_starts > 1) and, for the block bootstrap schemes, the out-of-bag
fraction -- each row gets a status ("ok"/"note"/"warning") and a
one-line suggestion.
evzinb() / evinb() gain family = , an evinf_family() object (or a
bare string, shorthand for evinf_family(count = string)) naming the
count-state distribution (count = ) and the zero process (zero = ).
count = "nbinom" (the default, reproducing today's model exactly) or
"poisson". A Poisson count state drops the dispersion parameter
Alpha.NB entirely -- not one fewer degree of freedom silently fixed at
a boundary, but genuinely absent from par.all, coef(), vcov(),
confint(), tidy(), and every bootstrap replicate's coef;
summary()/glance()/print() show it as absent rather than NA. The
Poisson M-step matches glm(family = poisson())'s closed-form MLE to
numerical precision (verified directly, and via the full
evzinb()/evinb() API).zero = "mixture" (the default) or "hurdle" (evzinb() only --
evinb() rejects it unconditionally, having no zero state to hurdle
over). Under a hurdle zero process the zero state owns every zero
outright (no mixing with the count state there, and the posterior for a
y = 0 row is forced to the zero state rather than derived from a
density ratio) and the count state is zero-truncated for y > 0 (its
log-density gains -log(1 - f0), f0 = P(Y=0) under the untruncated
count distribution; the new score/Hessian terms this requires were
derived symbolically and verified against central finite differences
before being committed). With the extreme-value state pinned
unreachable, a hurdle fit matches pscl::hurdle()'s coefficients and
log-likelihood to numerical precision, for both count = "nbinom" and
count = "poisson".residuals(), predict(type = "quantile"), simulate() and
classify_states() all dispatch on the fitted family (count
distribution and zero process alike); print()/glance() report it
(e.g. "poisson/hurdle").family = evinf_family()) fits are numerically unchanged
(audit §5.5, round9 E.0/E.1/E.2).compare_models() gains poisson_comparison / zip_comparison, mirroring
the existing nb_comparison / zinb_comparison: Poisson (glm()) and
zero-inflated Poisson (pscl::zeroinfl(dist = "poisson")) competitor
baselines. Both default to TRUE only when the fitted model itself used
family = "poisson", FALSE otherwise -- user-overridable either way;
zip_comparison is unavailable for evinb objects, with the same
defaults-to-FALSE-with-a-message / errors-if-TRUE behaviour as
zinb_comparison. predict.poissonboot() / predict.zipboot(),
tidy.poissonboot() / tidy.zipboot() and
glance.poissonboot() / glance.zipboot() are registered alongside the
existing nbboot / zinbboot methods (neither reports a dispersion
parameter: alpha is NA in glance()); compare_fit() and
oob_evaluation() need no changes, having already treated their model
list generically (round9 E.3).
evzinb() / evinb() / add_bootstraps() gain bootstrap_scheme = for
panel and time-series bootstrapping (round9 F): "iid" (plain row
resampling, the default when block is not given), "cluster" (resample
whole block units -- what block = has always done, and still the
default when block is given), "moving_block" or "stationary"
(block-resample each unit's own time series; Kunsch 1989 / Politis and
Romano 1994). The two block schemes need a new time = argument (a time
index, required to be strictly increasing within every block unit's rows
as they already appear in the data -- never sorted for you) and, optionally,
block_length = (NULL, the default, uses ceiling(T^(1/3)) per unit,
messaged once at the original fit, not on every bootstrap replicate).
block is optional for the two block schemes: the whole dataset is treated
as one unit when it's omitted. Every bootstrap replicate gains
$oob_fraction next to $boot_id -- with overlapping blocks the
out-of-bag set is both smaller and more temporally correlated than under
i.i.d. resampling, so out-of-bag error (oob_evaluation()) is optimistic
relative to genuine forecasting performance under those two schemes.
lr_test(bootstrap = ), compare_models() and oob_evaluation() needed no
changes: all three already just consume whatever boot_id the resampling
generator returns. Default (bootstrap_scheme = "iid" / "cluster" via
block) fits and bootstraps are numerically unchanged.
offset() is now supported in formula_zi and formula_evi (already
supported in formula_nb), not just the count component. An offset there
shifts the corresponding multinomial-logit log-odds against the count
state -- e.g. formula_evi = ~ x + offset(log(population)) for a
probability of an extreme event that scales with exposure. formula_pareto
still rejects offset(): a shape parameter has no clear reading for one.
There is no inheritance: a component that defaults to formula_nb gets
that formula with any offset() term stripped, exactly as before -- an
offset applies only where it is written explicitly. Offsets are carried
through the bootstrap (indexed by boot_id, like every other design
column) and through predict(newdata = ) (which requires the offset
variable, like the existing count-component offset already does). Default
(no-offset) fits are numerically unchanged (audit §5.6, round9 D.1).
evzinb() / evinb() gain weights = , given as a bare column name, a
string naming a column, or a numeric vector. Frequency-weight
semantics: every observation's contribution to the log-likelihood and to
the EM/M-step accumulations is multiplied by its weight (a fit with
integer weights matches a fit on the row-duplicated data to machine
precision at the likelihood level, and to 1e-6-1e-8 in the fitted
coefficients -- floating-point accumulation order can shift which exact
point the EM's own convergence tolerance stops at, the same phenomenon as
round8 A.3). nobs() -- and therefore AIC, BIC and the approximate
t-based p-values -- now use sum(weights), not the row count;
glance() gains sum_weights alongside the unchanged, row-count nobs,
so both are visible together. Weights must be positive and finite;
non-integer (analytic) weights are allowed, though the AIC/BIC
interpretation above assumes frequency weights. weights = does
not give design-based standard errors for survey data -- for that,
resample primary sampling units with block = instead. The bootstrap
resamples rows exactly as without weights and carries each drawn row's
weight along (audit §5.6, round9 D.2).
texreg support (round10 I.5, audit §5.9, decision D5): skipped.
modelsummary::modelsummary() (?gof_map_evinf, ?tidy.evzinb) already
covers regression tables for evzinb/evinb models; revisit only if a
user specifically asks for texreg.Implements section 4 of the internal package audit (new functionality) and the
review follow-ups in dev/review_round1.md.
$degenerate / $degenerate_reason): a replicate
whose EM did not converge; one with a non-finite fitted coefficient, or one
larger in absolute value than evinf_control(coef_limit = ) (default 50,
on the linear-predictor scale); or one whose smallest fitted Pareto shape is
non-finite or below evinf_control(alpha_floor = ) (default 0.001).
failed_bootstraps() lists them (new type column, "error" /
"degenerate"); glance() gains n_degenerate_bootstraps and n_bootstraps
now counts only usable replicates; print() / summary() report the count.
Every bootstrap-summary function (coef() bootstrap extractor, vcov(),
confint(), tidy(), summary(), predict(), marginal_effects())
excludes degenerate replicates by default and takes
exclude_degenerate = FALSE to keep them.C_EV estimate lands on an endpoint of the
candidate grid is \emph{not} treated as degenerate (dropping those would bias
the bootstrap distribution of C_EV inward); instead the count is stored as
object$n_c_on_boundary, reported by glance() and the print() boundary
note, and a single warning() suggests widening c.lim.evinf_control() bundles the EM tuning settings; evzinb() / evinb() gain
a control argument. The individual tuning arguments still work but are
deprecated in favour of control.c.lim = NULL); it is printed with a message. c_profile() /
plot_c_profile() show the log-likelihood profile over that range;
glance() gains n_above_c and n_em_steps; the fitted object carries
$c_profile, $c_trace, $loglik_trace.offset() in the count-component formula is supported
(mu_NB = exp(x'b + offset)), e.g. y ~ x + offset(log(exposure)).coef(), vcov() (bootstrap covariance), confint(),
logLik() / AIC() / BIC(), nobs(), formula(), terms(),
model.frame(), fitted(), residuals() (response and randomized quantile),
simulate(), update(). revzinb_fit() / revinb_fit() are superseded by
simulate().add_bootstraps() extends an existing fit with more replicates;
failed_bootstraps() returns the error messages of the replicates that
failed. Bootstrap seeds are recorded in $boot_seeds.lr_test(bootstrap = TRUE) and compare_models()) through progressr (shown when verbose = TRUE or a progressr handler is set), replacing the old cat() messages.classify_states() returns the prior and posterior state classification per
observation; state_table() cross-tabulates them.predict_grid() builds a one-variable prediction grid (other covariates held
at their mean/median or modal level, or pinned via at) and evaluates the
model over it in long form, optionally with bootstrap confidence intervals.plot() methods for fitted models: type = "states" (prior state
probabilities over a covariate), "prediction" (harmonic-mean prediction with
a bootstrap ribbon and quantile lines), "coefficients" (bootstrap
coefficient densities), "ppc" (observed-vs-expected binned frequencies) and
"ppc_quantiles" (observed vs. simulated tail quantiles with a 5-95% band).compare_fit() summarises the paired bootstrap differences (evinf - compared, negative favours the extreme-value model) in AIC, BIC and
out-of-bag RMSE / RMSLE for an evzinbcomp object, with the proportion of
bootstraps favouring the extreme-value model;
plot(), tidy() and glance() methods for evzinbcomp and an
oob_evaluation() method that tabulates the out-of-bag error per model.marginal_effects() computes average marginal effects (central difference for
numeric covariates, level-vs-reference contrasts for factors) on the harmonic
mean, the state probabilities or a predicted quantile, with bootstrap
confidence intervals. For type = "quantile", whose bootstrap distribution is
heavy-tailed, the reported std.error is a robust scale from the percentile
interval rather than the raw standard deviation.marginaleffects compatibility: evzinb / evinb models register with
marginaleffects on load, so marginaleffects::avg_slopes() and friends work
with bootstrap delta-method standard errors.gof_map_evinf() returns the modelsummary goodness-of-fit map (with an
extra argument to append rows); the bundled gm_evzinb data object is now
generated from it and gains an obs_above_c_ev row.evzinb() / evinb() accept block as a bare column name
(block = id) as well as a string.print(summary(model)) reports the number of observations at or above C_EV
and notes when the log-likelihood was recomputed.predict(type = "quantile") and marginal_effects(type = "quantile")
parametrised the negative-binomial part of the mixture by Alpha.NB where
every other part of the package (the likelihood, residuals(type = "quantile"), simulate()) uses size = 1 / Alpha.NB. Predicted quantiles
now use the correct dispersion; the effect is small for evzinb (Alpha.NB
near 1) and can be large for an over-dispersed evinb fit.
marginal_effects(type = "quantile") no longer returns zeros. The mixture
quantile is rounded to an integer for predict(), but a central finite
difference of that step function is zero almost everywhere; the marginal-
effects path now differentiates the continuous (unrounded) quantile.
quantiles_from_*() gain a round argument (TRUE by default;
predict(type = "quantile") is unchanged).
marginal_effects(method = "difference") is now implemented for numeric
covariates (average change from a delta-unit increase, default 1 unit); it
was previously accepted and ignored. type = "quantile" defaults to
method = "difference".
marginal_effects(type = "quantile") caps the rows it averages over at
n_max (new argument, default 500; automatic above nrow * n_bootstraps = 2e5, or whenever set explicitly; n_max = Inf disables it), with a message.
The per-observation mixture-quantile solve made this prohibitively slow on
data the size of hks.
update() no longer embeds the whole data frame in the refitted model's
stored call. evzinb() / evinb() keep data as the expression the user
passed; update() re-evaluates it, falling back to the embedded copy (bound
to a symbol) only when the expression can no longer be resolved. model$call
and str(model) stay small on large data.
evzinb() / evinb() resolve a bare block name to the matching column of
data before trying to evaluate it, so an unrelated object of the same name
in the calling environment no longer shadows the column.
add_bootstraps() errors if boot_seed was already used for this model (the
%dorng% stream is fully determined by the seed, so it would silently
duplicate existing draws and shrink the bootstrap SEs). evzinb() / evinb()
now always record the seed actually used (drawing one when boot_seed = NULL),
so object$boot_seeds is complete.
mixture_p() (the mixture CDF behind residuals(type = "quantile")) now
uses the discretised Pareto CDF that matches the pmf the likelihood uses
(1 - (C/(y+1))^a for y >= C), rather than the continuous 1 - (C/y)^a.
Randomized quantile residuals shift by a small amount in the extreme-value
tail (on genevzinb2, 12 of 100 residuals move, max change 0.015).
evinb() no longer runs the EVZINB EM step (which updates the zero-inflation
block) during the warm-up phase -- a copy-paste artefact from the pre-0.10.0
code. The zero-inflation intercept is now held fixed throughout, as the model
intends. Estimates are unchanged to machine precision on the tested data (the
fixed intercept starts far enough from zero that the spurious update was
numerically inert).
Block-bootstrap: the block variable is now included in the single na.omit()
step, so a model whose block column has missing values where the model
variables do not no longer mis-indexes the resamples (previously could error
or silently sample the wrong rows).
lr_test() refits the restricted models with the full model's control
settings (candidate range for C, tolerances, pdf.pl.type, ...) and warm-
started from the full-model estimates. Previously it used the package
defaults, so a model fitted with a non-default c.lim was compared against a
restricted fit that searched a different candidate set for C_EV (the LR
statistic could even come out negative).
The internal check that the reported log-likelihood equals the log-likelihood
at the returned parameters no longer warns from inside every bootstrap; it
recomputes silently, records object$loglik_recomputed, and (only for the
full-sample fit, only when verbose = TRUE) emits a single message.
A user-registered parallel backend (doParallel, doFuture, ...) is left
untouched when a function is called with multicore = FALSE.
A non-finite value produced by an in-formula transformation (e.g. log(x)
with x <= 0) now raises an informative error naming the column instead of
misbehaving in the C++ code.
User-supplied init.Beta.NB is used (regression-test hardened).
evinb parameter count and predict.*boot(pred = "original") indexing fixes
(carried over from the 0.9.4 audit follow-ups).
predict() on the zinb slot of a compare_models() result with
type = "counts" (or type = "all", or type = "counts", confint = TRUE)
no longer errors with "$ operator is invalid for atomic vectors".
lr_test(bootstrap = TRUE) no longer errors ("Tibble columns must have
compatible sizes") when any bootstrap replicate failed or was flagged
degenerate. It now filters to the usable replicates first, like every other
bootstrap summary, and gains exclude_degenerate = TRUE (audit0.10 §1.1);
the results also report n_bootstraps_used next to n_failed_bootstraps.
summary(standard_error = FALSE) / tidy(standard_error = FALSE) no longer
error ("Column 'se'/'std.error' not found") when the model has bootstraps:
approx_t_value is now silently forced to FALSE along with it, and
print.summary.*() renders correctly with the se/approx_t columns
missing (audit0.10 §1.2).
compare_models(razorize = TRUE) refits NB/ZINB bootstraps on the wrong
rows whenever the razorised data had fewer rows than the full data (which
it always does): the bootstrap resample indices were used as-is against the
smaller data_razor, mixing in unrelated or out-of-range rows. Resamples
are now mapped onto data_razor's own rows before refitting. Winsorised
out-of-bag error is now computed against the raw (unwinsorised) outcome,
not the winsorised one, so it is comparable to the EVZINB/EVINB OOB error.
inner_nb()/inner_zinb()/boot_refit_one()/oob_evaluation()'s internal
try()s no longer print to the console on a failed replicate (audit0.10
§1.5).
marginal_effects(variables = ) errors on an unrecognised covariate name
(listing the available ones) instead of silently returning a 0-row tibble.
A numeric covariate that only enters the model's formulas wrapped in
factor()/as.factor()/cut() (e.g. y ~ factor(g_num)) is now treated
as categorical (level-vs-reference contrasts), instead of being perturbed
as a number and producing unseen factor levels; one wrapped only in
poly()/ns()/bs() keeps the numeric derivative/difference path with
its perturbation clamped to the covariate's observed range. A covariate
used both ways across formula components errors clearly (audit0.10 §1.10).
A singular M-step Hessian no longer aborts the fit (or, on some BLAS
backends, silently returns a wildly ill-conditioned step): the Newton step
now uses arma::solve(..., solve_opts::no_approx), which reports failure
instead of throwing or returning garbage; on failure it retries once with a
small ridge, and if that also fails the affected block keeps its
pre-step value for that EM step (audit0.10 §1.3). This is also the usual
reason a bootstrap replicate failed.
The EM outer loop (the C_EV profile update) could in principle run forever
if the profile oscillated between two candidate values. evinf_control()
gains max.c.iter (default 50), capping each phase separately; hitting the
cap in the convergence phase sets converge = FALSE (see the new
$c_converged / glance() column to tell this apart from the inner EM not
converging) and, for a full-sample fit, a warning() names the last two
C_EV values visited. em_profile_c() also now uses which.max(), so an
exact tie in the profile no longer returns a length > 1 c_hat and every
candidate being NaN/infinite errors clearly instead of silently breaking
the outer loop (audit0.10 §1.4).
The likelihood could underflow to -Inf for a large count or a small
Pareto shape: the discretised Pareto log-pmf now uses a cancellation-free
form (a*log(C/y) + log(-expm1(a*log(y/(y+1)))) in place of
log((C/y)^a - (C/(y+1))^a)), the NB/Pareto mixture is combined with a
log-sum-exp instead of summing on the natural scale, the NB log-pmf uses
lgamma() instead of an O(y) loop, and the multinomial state probabilities
(in the C++ M-step and in prob_from_evzinb() / prob_from_evinb()) use a
numerically stable softmax. Coefficients move by at most ~2e-10 and
log-likelihood by at most ~7e-12 on the identity fixtures -- well under the
1e-6 the round's numerical-identity gate would have required flagging
(audit0.10 §1.11).
evzinb() / evinb() now validate the response after na.omit(): a
negative, non-finite or non-integer value errors, naming the count and the
problem, instead of being silently misread (a non-integer count is
effectively ceiling()'d by the NB part of the C++ code but used exactly by
the Pareto part; a negative one produced a cryptic "argument is of length
zero" several layers down). evinf_control(prune.c.range = ) now requires
[0, 1) (was [0, 1]; 1 passed validation but crashed with "invalid
'size' argument"); em_c_candidates() clamps its internal sample size to
>= 0, skips pruning entirely below 3 candidates, and no longer risks
sample()'s length-1-population trap (sample(2, ...) samples from 1:2,
not the single candidate 2) on a 3-value grid. Pruning's random draw is
now reproducible without touching the caller's .Random.seed, via a
generalised evinf_seeded_sample() (audit0.10 §1.9).
lr_test() no longer emits evzinb()/evinb()'s deprecation warning on
every restricted refit (once per model, or once per model x bootstrap with
bootstrap = TRUE). The restricted refits now pass a modified copy of the
full model's control object (its init.* fields set from the full-model
estimate) instead of ~20 individual tuning arguments, which also means a
setting added to evinf_control() after the model was fitted (e.g.
max.c.iter) is carried over automatically instead of needing to be added
to a hand-picked argument list (audit0.10 §1.7).
evinf_control(pdf.pl.type = "exact") was accepted and stored but had no
effect: the M-step always used the continuous-Pareto gradient/Hessian for
the Pareto block, giving bit-identical estimates to "approx". "exact"
now uses the analytic gradient/Hessian of the discretised Pareto log-pmf
(the pmf the likelihood itself always uses) in the Newton step, computed in
the same cancellation-free form as the likelihood (audit0.10 §1.11). The
default stays "approx", and default estimates are unchanged. See
?evinf_control for what each option does (audit0.10 §1.8).
A bootstrapped p-value of exactly 0 (no draw crossed the estimate) is now
floored at 1 / B, where B is the number of usable bootstrap replicates
-- the true p-value could be anywhere in [0, 1/B), so reporting exactly
0 overstates the precision B draws can deliver.
print.summary.evzinb() / print.summary.evinb() show this as "< 1/B"
(e.g. "<0.01" for 100 usable bootstraps) instead of the default,
misleadingly precise "<2e-16" (audit0.10 §1.11).
update(model, data = new_data) kept the old data's data-driven
control$c.lim / control$init.C (both were written into object$control
at the original fit), silently reusing a candidate range for \eqn{C_{EV}}
chosen for data the model no longer uses. If the original range was itself
data-driven (not pinned by the user via evinf_control(c.lim = )) and
data is among the arguments being changed, they are now reset to NULL
and re-resolved from the new data, with the usual message.
?update.evzinb also documents that changing formula_nb. does not
propagate to a component formula that was originally left NULL and
inherited from it -- that component keeps its own formula regardless
(audit0.10 §1.11).
object$props / object$resp (and the fitted$prob_* / posterior_*
vectors derived from them) came from the E-step at the start of the last
EM step, one step behind the returned coefficients. They are now
recomputed at the actual returned parameters and final C_EV before
evzinb() / evinb() return, using the same stable softmax as the C++
E-step. Point predictions (predict(type = "harmonic") etc.) are
unaffected -- only these diagnostic quantities change, by up to ~0.008 on
the bundled example data (audit0.10 §1.11).
plot() on a fitted model errored on ggplot2 < 3.5.0 (scale_*_continuous()
didn't gain the transform = argument until 3.5.0). The four log1p-scale
call sites in plot.evzinb() / plot.evinb() now pick transform = or the
older trans = based on the installed ggplot2 version, and the (>= 3.5.0)
floor on ggplot2 in Suggests is removed (round8 0.1).
pdf.pl.type = "exact"'s gradient/Hessian went to NaN for a sharply
peaked Pareto block (large alpha * |log(C/y)|, e.g. C = 100, y = 1e4,
alpha ~ 403) because pareto_exact_derivs_fun() formed (C/y)^alpha and
(C/(y+1))^alpha separately in its numerator even though the denominator
already used a cancellation-free form; both underflow to exactly 0 there,
giving 0/0. The NaN failed safe -- the M-step just stopped updating the
Pareto block silently -- but is now fixed by factoring the same term out of
the numerator (review §3, round8 0.2).
object$fitted$y.hat.pl_* (the state-probability-weighted point
predictions) were computed from the E-step's one-step-stale prior state
probabilities, while object$props / object$resp (and the
fitted$prob_* / posterior_* vectors derived from them) already used the
final, recomputed ones -- so fitted$y.hat.pl_E.inv.y and
predict(type = "harmonic") disagreed by a small but nonzero amount
(cor() = 0.9999999, not 1). Both are now computed from the same final
props (review §4, round8 0.3).
oob_evaluation() on an evzinbcomp object masked a degenerate evinf
bootstrap replicate's out-of-bag error to NA in the evinf column only,
leaving the compared nb / zinb columns unmasked at that row -- so a
column-wise na.rm = TRUE summary (as the vignette's model-evaluation
section does) compared medians computed over different replicate sets.
Every column is now masked at the union of NA positions, and the row
mask is exposed as attr(out, "excluded") (review §5, round8 0.4).
compare_fit()'s aic / bic rows for a *_razor or *_winsor slot are
not comparable to the evinf model's -- a razorised fit is estimated on
fewer observations, and a winsorised fit on a different outcome -- even
though the rmse / rmsle rows now are (OOB error is always against the
raw outcome). Those two metrics now come back NA for those slots, with a
one-line footnote from print() when any are present (review §6, round8
0.5).
The warm-up phase's C_EV profile could hit max.c.iter without settling,
same as the convergence phase, but nothing recorded it -- $c_converged
and its warning() only ever covered the convergence-phase loop. Recorded
now as $c_warmup_capped, with its own warning() for a full-sample fit;
it does not affect $converge (review §7, round8 0.6).
inner_nb() / inner_zinb() (the per-bootstrap refits behind
compare_models()) computed out-of-bag error against data[-boot_id, ],
which silently returns zero rows rather than all of them when boot_id is
empty -- a resample where razorising leaves no drawn row behind. Both now
guard explicitly and fail the same way a caught MASS::glm.nb() /
pscl::zeroinfl() error already does, so every existing
inherits(b, "try-error") check picks it up (review §7, round8 0.7).
CI: added an ubuntu-latest / oldrel-1 job to R-CMD-check.yaml -- the
job that would have caught round8 0.1's ggplot2-version failure (round7
prompt E.4, review §2, round8 0.8).
The count-likelihood closed form lgamma(y + 1/alpha) - lgamma(1/alpha) - lgamma(y + 1) (round 7) partially cancels for a large y with a small
1/alpha -- checked against a Kahan-summed reference loop (plain
summation of the pre-round-7 loop is not itself a trustworthy reference at
this scale) at up to ~1.5e-10 relative error for y up to 1e6 and alpha
down to 1e-3. Switched to the algebraically identical, cancellation-free
-log(y + 1/alpha) - lbeta(y + 1, 1/alpha), which was at least as
accurate in every case checked (review §7, round8 0.9). The Kahan-summed
reference used to check this also showed that the pre-round-7 loop
(plain summation) itself drifts by ~1e-8 at large counts -- so a fit of
large counts from before 0.10.0 had a slightly wrong log-likelihood all
along, not merely a slower one (round9 0.6).
A fitted Pareto shape (alpha_pl) that collapses toward 0 made
predict(type = "explog") return Inf and predict(type = "harmonic")
return an astronomically large finite number (e.g. ~1e22 x C on hks
with an ordinary specification) instead of erroring or warning -- both
involve exp(1 / alpha_pl) or 1 / alpha_pl. A shared helper
(evinf_clamp_alpha_pl()) now floors alpha_pl at
evinf_control(alpha_pl_floor = ) (default 0.01) in harmonic_calc(),
explog_calc(), the derived Pareto-tail summaries in $fitted
(em_fitted_values()), and the quantile-prediction clamp that already
existed for quantiles_from_evzinb() (extended to quantiles_from_evinb(),
which previously lacked it). The helper warns once per call, naming how
many observations were clamped -- except inside em_fitted_values(), which
runs on every EM fit including every bootstrap replicate, so it clamps
silently there rather than warning on routine fits. glance() gains
min_alpha_pl (the true, unclamped value) and print() adds a note when it
falls below the floor, so a collapsed extreme-value shape stays visible
without digging (review §2, round9 0.1). The opposite direction --
alpha_pl overflowing exp() to literal Inf (an extreme Beta.PL
coefficient on a sparse factor level, rather than a collapsing tail) --
isn't a case a floor helps with: harmonic_calc()'s
(1 + alpha_pl) / alpha_pl computed Inf/Inf = NaN there. Rewritten as
1/alpha_pl + 1, exact for any finite alpha_pl and correct at
alpha_pl = Inf too (the extreme-value contribution's mathematical limit
there is C), so no ceiling is needed either.
The round8 A.2 closed form for the NB dispersion Hessian's alpha-alpha
entry switched to its O(y) loop fallback only for y <= 30, but the
cancellation it guards against is governed by the product alpha_nb * y,
not y alone -- the near-Poisson regime (small alpha_nb) is exactly
where this mattered and exactly where the round8 grid (alpha >= 1e-3)
never checked. A fresh grid down to alpha = 1e-6 found the closed form
18% off a loop reference at y = 31, alpha = 1e-6 (alpha_nb * y = 3.1e-5), 0.47% off at y = 50, and still 0.27% off at y = 100 --
errors a y > 30 guard let straight through. The guard now switches on
alpha_nb * y > 0.1, which keeps a comfortable margin above where the
closed form's relative error crosses 1e-9 (observed worst case ~1e-11
once alpha_nb * y clears that cutoff); a large y together with a
vanishingly small alpha_nb still takes the O(y) loop rather than the
(still available, but not implemented here) small-product series
expansion -- accepted per the round9 brief, since the loop is only slow,
not wrong. This matters directly for the round9 0.11.0 Poisson count
family (alpha_nb -> 0 is its limit) (review §4, round9 0.4).
glance()$n_bootstraps no longer counts every bootstrap replicate --- it
counts the usable ones (neither errored nor degenerate). The three
columns n_bootstraps, n_failed_bootstraps and n_degenerate_bootstraps
now partition the number of replicates requested. The gof_map_evinf() row
labels were retitled accordingly ("Usable bootstraps", "Failed bootstraps",
"Degenerate bootstraps") and given title case throughout; the raw names
(i.e. the glance() column names) are unchanged.
The parallel backend is now future / furrr instead of
foreach / doParallel / doRNG (which are no longer imported).
Set a plan once for your session — future::plan(future::multisession, workers = 8) — and every bootstrap loop (evzinb(), evinb(),
add_bootstraps(), lr_test(bootstrap = TRUE), compare_models(),
predict(), marginal_effects()) uses it. The multicore argument now
defaults to NULL (respect the current plan); multicore = TRUE still
works as a shortcut that sets a temporary multisession plan for the one
call and restores the previous plan on exit, and multicore = FALSE still
forces sequential execution. See the new "Parallel processing" section in
?evzinb.
Because the RNG is now derived with furrr's per-element L'Ecuyer streams,
the bootstrap resamples drawn for a given boot_seed differ from earlier
versions (the resampling distribution is unchanged, and results no longer
depend on the number of workers). Pin results you need to reproduce.
future and furrr are new hard dependencies; foreach, doParallel,
doRNG and mistr are dropped.
The likelihood, CDF, quantile prediction, residuals, simulate() and the
posterior state probabilities now all share a single discretised Pareto
(R/dist_pareto.R): pmf \eqn{(C/y)^\alpha - (C/(y+1))^\alpha} for integer
\eqn{y \ge C}. The 'harmonic'/'explog' point predictions still use the
continuous-Pareto harmonic/geometric-mean formulas as a closed-form
approximation (documented in ?predict.evzinb), so they are the one place
that is not computed from the discretised distribution. Previously
predict(type = "quantile") / quantiles_from_*() /
marginal_effects(type = "quantile") built the mixture with mistr and
inverted it there. The new path bisects the same mixture CDF used by
residuals(type = "quantile") and the likelihood, so it also corrects the
negative-binomial dispersion those quantiles used (they parametrised the NB
by Alpha.NB rather than the canonical size = 1 / Alpha.NB). For an
evzinb fit, where Alpha.NB is near 1, predicted quantiles move by at most
a count or two; for an evinb fit with a strongly over-dispersed count
component the old quantiles could be off by a large factor. simulate() /
revzinb_fit() / revinb_fit() draw the extreme-value part with
floor(C * U^(-1/alpha)) (was round(...)), so simulated extreme-value
counts are ~0.5 lower on average and never below C.
The minimum R version is now 4.1.0. The NAMESPACE uses delayed S3 method
registration for insight / marginaleffects, which needs R >= 3.6.0; 4.1
is the common floor for packages using that pattern.
The approximate bootstrap runtime estimate is printed only when
verbose = TRUE (it was always printed before 0.9.4; this note was missing
from the 0.9.4 changelog).
predict(type = "all") and every confint = TRUE output now lead with the
canonical state-probability columns (pr_zero, pr_count, pr_evi), keeping
the deprecated pr_zc / pr_pareto as trailing duplicates — matching what
predict(type = "states") already did.
Default fits change slightly: with c.lim = NULL the candidate set for C_EV
is data-driven and init.C defaults to the median of that set, so estimates
from a default 0.9.x call are not bit-identical to a default 0.10.0 call. Pin
evinf_control(c.lim = ..., init.C = ...) to reproduce an older fit.
gm_evzinb gains a row (obs_above_c_ev) and its parameter row is now
formatted with 0 decimals; regenerate any cached copy with gof_map_evinf().
The hks column brv_AllLag is renamed brv_AllLag_log (it is on the log
scale, like the other _log columns). There is no alias — code referring to
the old name will get an immediate "object not found" / "column not found"
error.
predict() on the nb / zinb slot of a compare_models() result no longer
accepts type = "evinf" (a plain NB / ZINB has no extreme-value state);
match.arg() now rejects it.
The three internal marginal.effect.* helpers (never exported, superseded by
marginal_effects()) were removed.
compare_fit() sign convention. The paired bootstrap difference is now
evinf - compared (was compared - evinf); prop_evinf_better is now
mean(diffs < 0). The print header and the plot() x-axis label already
said "negative favours evinf" — it was the computation that disagreed with
them (a probe on genevzinb2 showed a positive median with
prop_evinf_better = 1 for a model that clearly wins on AIC). compare_fit()
/ plot.evzinbcomp() also gain exclude_degenerate = TRUE, and
oob_evaluation() now honours it for a single evinf model.
The default conf_level for predict() and marginal_effects() is now
0.95 (was 0.9), matching confint() and tidy(). Pass
conf_level = 0.9 explicitly to keep the old default.
Deprecations planned for removal in 0.11.0: the pre-0.9.4
component = "nb"/"zi"/"evinf" aliases (use "count"/"zero"/"evi");
passing the individual EM tuning arguments (max.diff.par, c.lim,
init.C, ...) directly to evzinb()/evinb() instead of through
control = evinf_control(...); and the deprecated pr_zc/pr_pareto
column names and fitted$prob_pareto/fitted$posterior_pareto duplicate
fields (use the canonical pr_zero/pr_evi and
prob_evi/posterior_evi). All of these still work in 0.10.0 and most
already warn when used. revzinb_fit()/revinb_fit() are superseded by
simulate() but are not deprecated and have no planned removal.
?vcov.evzinb, ?confint.evzinb and the new ?marginaleffects-methods page
note that c_ev is estimated on a discrete grid, so continuous delta-method
standard errors that perturb it should be read with care;
marginal_effects() (bootstrap) is the recommended route for effect
uncertainty. ?confint.evzinb gains an example using c_ev / alpha_nb.pkgdown reference index grouping the exported functions (Fitting,
Summaries & tables, Prediction & effects, Diagnostics & plots, Model
comparison, Simulation, Data).?hks now documents the columns as they actually ship (the previous help
page listed several columns that are not in the data and omitted the _log
transforms) and explains that log1p(troopLag) etc. can be used directly in
the formula since 0.9.4.R CMD check, test coverage and the pkgdown
site; coverage and check badges in the README.compare_models() are replaced by a single boot_refit_family() helper
(output unchanged).R/a019.R is split into six documented,
de-duplicated files (R/em_candidates.R, R/em_profile_c.R, R/em_step.R,
R/em_fixed_c.R, R/em_fitted.R, R/em_fit.R), one EM/ECME driver
em_fit() for both models instead of the two ~90%-duplicated
*.regression.fun() pairs. Estimates are unchanged (a numerical-identity
test suite pins the pre-split output). The three unused
prediction.*.fun() legacy helpers were removed.evinf_control(prune.c.range = ...) now also thins the \eqn{C_{EV}}
candidate set for evinb() (it was silently ignored there before), and the
"more than 100 candidates" warning now fires for evinb() too.?hks: the _log columns are documented as a rule (log1p() for the
personnel counts, log() for the episode duration), and the help notes that
the bundled data carry no conflict identifier.?marginal_effects / ?tidy.evzinb / ?gof_map_evinf gain guidance on
harmonic-mean effect intervals and on the modelsummary shape argument.em_profile_c() used to call
log_lik_fun() once per C_EV candidate, redoing the state-probability
and count-log-likelihood computation (which doesn't depend on C) every
time. New C++ log_lik_profile_fun() computes that shared work once and
profiles the whole candidate grid in one call; em_profile_c()'s return
shape, which.max() behaviour and all-infinite error are unchanged.
Verified to 1e-10 against the old per-candidate calls, and the selected
c_hat is unchanged (not just close) on every identity fixture.delldtheta_nb_i_fun() / d2elldtheta2_nb_i_fun() (the NB dispersion
score and Hessian) are closed forms via digamma()/trigamma(). The
Hessian's closed form cancels catastrophically for a small alpha * y
(not just a small y), so it's used only above a cutoff on that product,
with the loop (already trivially cheap in that regime) kept below it --
see round9 0.4 below, which corrected the cutoff variable this bullet
originally (and wrongly) described as y > 30. Verified to 1e-9 against
the old loops over y in {0, 1, 5, 100, 1e4, 1e5} x alpha in
{1e-3, 0.01, 0.5, 1, 5, 50}; the M-step gradient/Hessian are unchanged
on the identity fixtures.log_lik_fun() and
update_bfgs_fun()'s per-observation loops rebuilt a linear predictor
(trans(X.submat(i, ...))) and, for the four M-step BFGS blocks,
allocated a fresh gradient/Hessian block via a per-row function call --
on every row, every BFGS sub-iteration, and log_lik_fun() alone is
called dozens of times per em_step() (once per evaluation of each of
the four stats::optimise() line searches in R/em_step.R). Replaced
with precomputed linear predictors (one matrix-vector product per block)
and gradient/Hessian accumulation as X' w / X' diag(w) X. A full
hks fit went from 1.1x faster (A.1 + A.2 alone) to a platform-dependent
1.9x-8.4x faster (round9 0.5, review §3: the original single 8.4x/163x
headline was measured on Apple Accelerate and did not reproduce on
Linux/OpenBLAS, where the same hks fit was only 1.9x faster -- quote a
range, not one number, and see inst/bench/bench_evinf.R, now tracked and
CI-runnable, for how to reproduce either end of it); see A.4's benchmark
table for the Apple Accelerate figures. Verified against the pre-A.3 code on 5 cases
(genevzinb2 and hks, both Pareto types, plus an evinb fit) at
machine precision (max abs diff ~2e-12); the identity fixtures are
unchanged at 1e-8, with c_hat exactly (not just closely) unchanged on
those cases. That comparison is on a single call at identical fixed inputs,
not on a full fit's trajectory: the reordering perturbs each Newton step at
the ~1e-12 level, and on data with a close second optimum that is enough
for the iterated EM to converge somewhere else. On one factor-covariate
fit this moved the log-likelihood from -253.5553 to -253.6816 and
min(alpha_pl) from 1.91 to 0.973 on some platforms (review §1) --
"estimates unchanged" above is therefore too strong as a blanket claim;
read it as "unchanged on every case checked so far", with a dedicated
fixed-trajectory regression test now guarding that specific case (round9
0.3) and the alpha_pl -> 0 consequence guarded separately (round9 0.1).This release works through sections 1-3 of the internal package audit: it fixes
substantive statistical issues, a number of usability bugs, and brings the
documentation and metadata up to date. A testthat suite has been added.
lr_test() now reports the standard
2 * (logLik_full - logLik_restricted) statistic (and uses it for the
bootstrap reference distribution). Statistics reported by earlier versions were
half the correct value; p-values change accordingly but the qualitative
conclusions in Randahl & Vegelius (2024) are unaffected. Degrees of freedom are
now the number of design-matrix columns removed across components (so removing
a four-level factor contributes three degrees of freedom)."count", "zero",
"evi" and "pareto" everywhere (tidy(), coefficient_extractor(),
summary()). The old names "nb", "zi" and "evinf" still work but emit a
deprecation warning. predict(type = ) additionally accepts "zero" for
"zi" and "evi" for "evinf".tidy() default is now component = "all" (was "zi"), so a bare
tidy(model) / modelsummary(model) shows every component. The y.level
values from tidy(component = "all") are zero, evi, count, pareto.predict(type = "states") columns are pr_zero / pr_count / pr_evi;
model$fitted$prob_pareto / posterior_pareto become prob_evi /
posterior_evi; the props / resp matrices are labelled zero / count /
evi.compare_models() returns the fitted model in $model; $evzinb remains
as an alias. It now also accepts evinb objects (with zinb_comparison
forced to FALSE).generics moved from Depends to Imports; tidy() and glance() are
re-exported, so library(broom) / library(generics) is no longer required.MLmetrics, stringr and stringi removed.compare_models() fitted the bootstrapped ZINB with the zero-inflation formula
for both model parts; it now uses the full count | zero formula, matching
the full-sample fit.summary() and tidy() no longer error on models fitted with
bootstrap = FALSE; they return the point estimates with a message.summary(p_value = "both") now returns both the bootstrap and the approximate
p-values (the approximate branch was previously unreachable).lr_test() works on evinb objects (it passed the wrong data object before).model.matrix() rather than
model.frame()[, -1], so factor and character covariates, interactions,
poly() / splines and in-formula transformations such as log(x) now work and
are labelled correctly. na.omit() is applied once to the union of all
component variables, keeping the components row-consistent.init.Beta.NB is no longer silently discarded when the NB and
EVI formulas differ in length.evinb() bootstraps in the returned object are now named.predict.nbboot() / predict.zinbboot() are registered as S3 methods.summary() gains a print method; print.evzinb() / print.evinb() show the
correct component formulas plus C_EV and the count of observations at or
above C_EV.glance() methods return unrounded values and gain converged,
n_bootstraps and n_failed_bootstraps columns; gm_evzinb updated to match.by" messages are emitted.run_evinb()'s default c.lim matches evinb(); pdf.pl.type is validated.logLik / AIC /
BIC (with a warning) if not.tidy(confint = ) is implemented (percentile or approximate intervals).revzinb_fit() / revinb_fit() draw Pareto values in a single vectorised call.evinb objects no longer count the (fixed) zero-inflation coefficients as free
parameters: glance()$npar, AIC, BIC and the approximate-test degrees of
freedom are corrected accordingly.predict() on nbboot / zinbboot objects with pred = "original" used a
stray loop index and errored; it now predicts from the full-sample model.inst/CITATION and inst/REFERENCES.bib updated to
International Studies Quarterly 68(4): sqae137 (2024),
doi:10.1093/isq/sqae137.hks: policeLag and militaryobserversLag are now correctly described as UN
police and military observers; all personnel counts are documented as being in
thousands.genevzinb / genevzinb2: "independent" (not "dependent") covariates;
"EVZINB" spelling.prune.c.range documented consistently across all four fitting functions;
c.lim reworded.man/ regenerated with roxygen2 7.3.3; src/a019_old.cpp renamed to
src/evinf.cpp; empty src/coef_extractor.cpp removed.\donttest{} and reduced to
n_bootstraps = 5, multicore = FALSE.README filled in.predict() can return the raw bootstrap draws (return_bootstraps = TRUE).prune.c.range for thinning the candidate set of C when it contains many
unique values.evzinb to evinf.revzinb_fit() / revinb_fit() for drawing from a fitted model.predict() changed.