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.