First release. iop estimates ordered probit and ordered logit models whose
outcome has an inflated category, together with the plain ordered models, on
one engine with one set of methods.
oprobit(), ologit(): ordered probit / logit, with survey weights,
offsets, partial proportional-odds (non-parallel) effects via parallel =,
and analytic, robust, or cluster-robust standard errors.iop(): inflated ordered probit for any single inflated category
(inflate = "bottom" | "middle" | "top" or a category label), the
zero-/middle-/top-inflated ordered probit of Harris and Zhao (2007),
Bagozzi and Mukherjee (2012), and Bagozzi, Joo and Mukherjee (2024), with
optional correlated errors (correlated = TRUE).iol(): the inflated ordered logit counterpart.y ~ x1 + x2 | z1 + z2, a C++
likelihood with analytic gradients (a deterministic Gauss--Legendre
bivariate-normal algorithm for the correlated form), column-scaled
quasi-Newton optimization with an exact-Hessian Newton polish, and one S3
class ("iord") for every method.vuong() and the Vuong block of inflation_test() are now consistent
with the package's frequency-weight semantics: the statistic equals the
one computed on duplicated rows (weight total as the sample size, weighted
mean and variance of the per-row log-likelihood difference); previously
the unweighted row count and the unweighted standard deviation of the
weighted contributions were used, which gave a different statistic under
non-constant weights. The BIC (in summary(), glance(), texreg,
inflation_test()) uses the weight total as n for the same reason.classification()'s confusion table is weighted (sums of weights per
cell), so it agrees with the weighted share correctly classified printed
under it.w / P to Inf (and Inf * 0 to NaN) in the gradient: the divisor is
floored at the smallest normal double. The case arose with an extreme
offset and, in random-intercept fits with a large unit variance, at the
outer quadrature nodes of a high nAGQ, where it could leave the fit at a
non-stationary point with a warning.$max_grad is reported in the optimizer's standardized parameterization,
the units the convergence rule uses.print() and summary() state the response
levels in the order used.predict(newdata = ) returns NA rows, as documented, when an all-NA
covariate column arrives as logical (data.frame(x = NA)); previously that
case failed with "subscript out of bounds".?ame points to first_difference(newdata = ); inflation_test(boot = )
reports the expected running time when it is long and documents the seed
guarantee, the cost, and the independence of its simulated data from any
cluster structure; se = "bootstrap" documents the seed guarantee;
?update.iord notes what a fit restored with readRDS() needs;
?classification cautions about modal accuracy for minority categories;
the data pages state the provenance and terms of the bundled data sets.weights are frequency weights and that survey
(probability) weights call for se = "robust" or cluster =;
inflation_test(boot = ) notes that its parametric bootstrap draws
independent observations and does not reproduce a cluster structure.vignette("iop") gains two worked examples outside political science: a
zero-inflated ordered probit for the number of extramarital affairs
(AER::Affairs) and a middle-inflated ordered probit for self-rated health
in the 1988 NMES (AER::NMES1988); AER is in Suggests and the chunks run
only when it is installed.$fe_units lists the unit levels actually estimated (a unit that
disappears through missing covariates is no longer counted).iord_probs_cpp, the
random-intercept entry points) and for the two examples.fe =
variable's name (unitprice next to fe = "unit") is a covariate -- it is
printed, screened by the separation diagnostics, and reported by texreg
like any other; previously it was folded into the hidden dummies
($fe_index, print(), summary()$block, extract()).glance() gains converged, boundary, and ill_conditioned, and the
texreg GOF block gains "Converged" and "Boundary/ill-conditioned" rows,
so tables built from a saved fit carry the fit-quality flags.NA.fe_inflation = TRUE) are exempt from
the split-equation quasi-separation signatures, as the outcome-equation
dummies already were; a short panel with many inflation dummies no longer
gets a "simplify the inflation equation" message when the cause is the
dummies (the incidental-parameters warning, the conditioning flag, and the
NA standard errors already say so).split_test()'s score and
information (both forms) are now all computed in the optimizer's scaled
parameterization (every covariate column standardized) and mapped back,
so none of them depends on the units of the covariates. Previously the
Hessian was differenced in raw parameter units: a covariate in cents next
to a dummy could return a negative delta-method variance (an NA standard
error) for a perfectly identified fit, and split_test(information = "hessian") could report "not positive definite" on a well-behaved null
fit with a split covariate in dollars. Regression tests fit the same model
at two covariate scales.update() refuses a formula. with more than one |, and a two-part
formula. on a fit without an inflation equation, with a clear message;
the estimators refuse formulas with more than one |.split_test() documentation
(inst/mc/split_test_size.R, with its per-replication results; the log
is kept in the repository) ships with the package; the rank-deficiency message names the dependent column and says
so; the quasi-separation message says the estimate "may not be finite".ame() documents that it has no profile mode by design.ame() and first_difference(average = TRUE) now honour the fit's
weights when averaging over the estimation data (they equal the
duplicated-rows fit exactly); previously the averages were unweighted.$boundary also flags a
quasi-separated split (standardized split slope above 10 or intercept above
40, or degenerate split standard errors), not only the all-at-0/1 case;
$boundary_messages, $ill_conditioned (information matrix with reciprocal
condition number below 1e-12 in the optimizer's scaled parameterization, so
that covariate units play no role, or a generalized inverse needed), and $se_na
(non-positive delta-method variances, reported as NA rather than 0) are
stored and printed by summary() and print(); the thresholds are listed in
?iop ("Diagnostics and fixed thresholds").fe = is now explicit.split_test(information = "opg" | "hessian"): the outer-product LM of
Brown, Harris and Spencer stays the default (correctly sized under the null
in simulation, always defined); the observed-information form is available
and is reported as NA with a note when the Hessian is not positive
definite at the restricted estimates. The documentation now states plainly
that the LM is a local screening statistic that can be far from the LR
under a strong alternative (267 vs 57 on bp) and that the LR is the one to
report.update() method for "iord" fits: two-part formula. arguments
(. ~ . + x3 | ., . ~ . | . + z2) update the outcome and inflation
equations separately; a one-part formula. changes the outcome equation
and keeps the inflation equation (stats::update.formula() alone folds
the | into the outcome part).print(inflation_test(...)) notes when the parametric-bootstrap LR was not
run; vuong()/lr_test() check that the two fits use the same data rows,
not only the same response vector; update() semantics, missing-value
handling, and newdata rules are documented on the estimator pages; the
jackknife covariance is justified (Dhaene and Jochmans 2015) and noted in
summary(); maxit documentation clarified.$boundary now flags it (ordered-stage share
of the inflated category's probability mass below 1e-6), names the
cutpoint, and summary()/confint() report its standard error as NA
instead of the generalized-inverse value. $flat_hessian is now computed
from the final Hessian of every fit, not only when the Newton decrement
was unavailable.0 + g, an
interaction x:g without the x main effect next to poly(x, 2)) is
rejected with a message naming the columns, where previously it fit with a
singular information matrix; an explicit 0 + / - 1 in the outcome
formula is ignored with a message (factors are coded as with an intercept).pbvn_cpp() and friends guard empty inputs;
mundlak() no longer uses <<-; parallel listed in Suggests; the CI
workflow runs the full suite (NOT_CRAN = true).split = "category" in iop() / iol(): every non-inflated category gets
its own split (inflation) equation -- and its own rho with
correlated = TRUE -- the generalised GZiOP / GMiOP of Brown, Harris and
Spencer (2020), which nests the common-split ZiOP / MiOP / TiOP. The fit
starts from its common-split counterpart ($loglik_common); coefficients
are infl_<term>:<category> and rho:<category>; regime quantities become
one column per non-inflated category. riop() draws from the generalised
process when gamma is a matrix.split_test(): Lagrange-multiplier (score) test of the common split from
the common-split fit alone, plus the likelihood-ratio test from the
category-specific refit.inflation_test(boot = R): a parametric-bootstrap likelihood-ratio test of
the inflated against the plain ordered model (simulated from the fitted
ordered model), the boundary-aware alternative to the Vuong test; the help
page documents the Wilson (2015) / Dale and Sirchenko (2021) objection to
the Vuong test for this nested comparison.predict(type = "zeros"): the two components of the inflated-category
probability (through the inflation process / through the ordered stage,
Harris and Zhao's two types of zeros), with se.fit; first_difference()
and ame() gain decompose = TRUE to report effects on both components.
Also predict(type = "mean") (expected category index) and
type = "cumulative".classification(): confusion table, share correctly classified, Brier and
ranked probability scores, mean log score, and per-category precision,
recall, and adjusted noise-to-signal ratio (Dale and Sirchenko's
ziopclassification), on the estimation data or on new data.ame(stage = "outcome" | "inflation"): stage-specific average marginal
effects, matching first_difference(stage =).diord(), piord(), qiord(), riord(): the d/p/q/r family of the
(inflated) ordered response -- probability mass, cumulative probability,
quantile, random draws -- from explicit parameters (eta, tau, a, k,
link, rho, including the category-specific split) or from a fitted
model at newdata.se = "bootstrap" (with nboot, cores) on every estimator:
nonparametric bootstrap over rows, clusters (cluster =), or
random-intercept units, with relabelled duplicates; replicate estimates in
$boot; confint(type = "percentile"). Recommended by Dale and Sirchenko
for the error correlation in small samples.summary(), confint(), vcov(scale = "natural" | "internal");
se = "analytic" | "robust" | "cluster".vuong() (raw, AIC-, BIC-corrected), lr_test() for nested pairs,
inflation_test() (inflated vs plain ordered, refit internally),
compare_models().ame(): average marginal effects (derivatives for continuous covariates,
discrete changes for 0/1 covariates and factor levels) on every category
and on the regime probability, with delta-method intervals; plot()
methods for ame() and first_difference() tables.
predict(se.fit = TRUE): delta-method standard errors for predicted
probabilities.
parallel_test(): likelihood-ratio tests of the parallel-regression
assumption term by term and jointly; parallel = "auto" relaxes terms
automatically (forward autofit).
fe_correction = "jackknife": split-panel jackknife bias correction of
the common parameters in fixed-effects fits (Dhaene and Jochmans 2015).
predict() types: category probabilities, modal class, ordered-stage
probabilities, regime and inflation probabilities, posterior probability of
being an inflated case, linear predictors.
first_difference() by category and regime, by stage ("both",
"outcome", "inflation"), at a profile or averaged over the data, with
delta-method or simulation intervals.
simulate() for DHARMa-style simulated-residual diagnostics; riop()
draws data from the model's DGP.
Panels: unit random intercepts by adaptive Gauss--Hermite quadrature
(re =, optionally also in the inflation equation with re_inflation = TRUE; ranef() for the empirical-Bayes unit effects, marginal and
conditional predictions), unit fixed effects as dummies (fe =,
fe_inflation = TRUE; units without identifying variation dropped with a
message; incidental-parameters warning for short panels), and the
mundlak() correlated-random-effects helper (both equations). A Monte
Carlo (inst/mc/) documents the trade-off: with a covariate correlated with
the unit effect, unit dummies are biased +26/+11/+5/+2.5 percent at T =
4/8/16/32, a random intercept +21/+11/+6/+3 percent, pooled +12 percent,
Mundlak within 0.3 percent.
broom tidy()/glance()/augment(); texreg extract().
bp (Besley and Persson 2009 political-violence panel, as in Bagozzi et
al. 2015), pta (Baccini, Dur and Elsig 2015 escape-flexibility data, as in
Bagozzi, Joo and Mukherjee 2024), repression (Girod, Stewart and Walters
2018 repression-of-campaigns data, as in the same article's appendix).vignette("iop") (getting started), vignette("quantities")
(predicted probabilities, first differences, average marginal effects,
tables, diagnostics), vignette("panels") (random intercepts, fixed effects
and the jackknife, the Mundlak device, the Monte Carlo), and
vignette("model") (likelihood, identification, estimation, boundary cases,
tests, validation).?"iop-package") indexing every function; examples
on every method page; @family cross-references between related pages; a
pkgdown site configuration (_pkgdown.yml) and GitHub Actions workflows for
R CMD check and the site.With fe =, a covariate that does not vary within units is refused with a
message that names it (it is an exact linear combination of the unit
dummies); previously the generic rank check blamed an arbitrary dummy
column. The same check applies to the inflation equation with
fe_inflation = TRUE.
The boundary flag distinguishes its two inflation cases in the warning: a split that collapses to the plain ordered model (every regime probability at 1) and a quasi-separated split (a few units at regime probability 0, the rest at 1, coefficients not finite).
fe_correction = "jackknife" stores the two half-panel coefficient vectors,
log-likelihoods, sizes, and dropped-unit counts in object$jackknife, so
the correction can be inspected before it is reported.
Formulas with transformed terms (poly(), ns(), scale(), I(),
interactions, expression responses) are built from per-equation model frames
so that predict() on new data uses the fitted bases (predvars); constant
or collinear covariates are rejected with a clear message; unused factor
levels are dropped.
Inflated fits are multi-start (binary-model starts at several slope scales
and signs, flat high-regime starts, and the ordered baseline); the
log-likelihood reached from every start is stored in $start_logliks. This
finds the second likelihood mode of the Besley--Persson ZiOP/ZiOPC with the
full inflation equation (published ZiOPC2: -1,370, rho -0.912).
The degenerate-inflation diagnostic is based on the fitted regime probabilities (all at 0/1), not on coefficient magnitudes.
lr_test() refuses pairs that differ in link, inflated category, or
random-intercept structure.
pbivnorm and mvtnorm to
1e-13; analytic gradients agree with numerical differentiation in every
model configuration.oprobit()/ologit() reproduce MASS::polr, ordinal::clm, and
VGAM::vglm(cumulative(parallel = FALSE)) (coefficients, cutpoints,
log-likelihoods, robust and cluster-robust SEs via sandwich); their
random-intercept versions reproduce ordinal::clmm(nAGQ = 15).iop() reproduces Table 1 of Bagozzi, Hill, Moore and Mukherjee (2015) on
the bundled bp data (OP, ZiOP, ZiOPC, ZiOPC2 log-likelihoods and rho),
Table 1 and the marginal-effects figures of Bagozzi and Mukherjee (2012;
MiOP and MiOPC, to the third decimal), and the TiOP/TiOPC log-likelihoods of
Bagozzi, Joo and Mukherjee (2024; Table 1 and appendix Table A.5) -- the
original likelihood code of each paper evaluated at iop's estimates equals
iop's log-likelihood to machine precision -- and matches the Python
idcempy implementation (ZiOP, MiOP) where that implementation converges.