iop 0.1.0

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.

Estimators

• 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.
• All four share the formula interface 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.

Second referee panel (2026-08-22)

• 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.
• A cell probability that underflows to a subnormal double no longer sends 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.
• A character response is refused with a message naming the alphabetical order it would have received; 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”.
• A data frame with no complete cases gives a message naming the cause; ?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.
• Documentation states that 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.

Final pre-submission additions (2026-08-22)

• 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).
• Tests for the remaining low-level C++ guards (iord_probs_cpp, the random-intercept entry points) and for the two examples.

Referee round 3 (2026-08-22)

• Unit-dummy columns are identified from the design’s term assignment, not by a name-prefix match: a covariate whose name starts with the 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.
• Column scales treat a covariate as constant only when its SD is negligible relative to its magnitude (1e-10), so covariates measured in billionths are standardized like any other (the absolute 1e-8 floor left them unscaled); the column scales are computed once per fit.
• Documentation notes why degenerate split/outcome standard errors stay visible while a cutpoint at +/-Inf is reported NA.
• The inflation-equation unit dummies (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).

Referee round 2 (2026-08-22)

• Standard errors, the condition number, and 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.
• The hurdle-type boundary flag (ordered stage empties the inflated category) also fires when the ordered-stage share is below 1e-2 and the adjacent cutpoint’s standard error already exceeds 50, so a fit the optimizer stopped just short of the boundary is named rather than left to the generic large-standard-error note; the degenerate-split-SE and large-SE notes are on the standardized scale.
• The quasi-separation signatures (standardized coefficient above 10, or standardized standard error above 50) are also applied to the outcome equation of every model (unit dummies exempt).
• 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 |.
• The size simulation behind the 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”.
• Tests for the low-level C++ guards (wrong-length parameter vector); ame() documents that it has no profile mode by design.

Referee round 1 (2026-08-22)

• 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 diagnostics are graded and surfaced: $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”).
• The Newton polish and the covariance use a central-difference Jacobian of the analytic gradient (2p gradient evaluations) instead of numDeriv’s Richardson extrapolation (about 8p), a 3-4x saving that matters for fits with many unit dummies; the documented scaling of 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.
• A boundary case the split-side diagnostics did not name: the ordered stage can empty the inflated category (every such observation attributed to the inflation process, a hurdle-type split), leaving the adjacent cutpoint at +/-Inf with no information; $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.
• The outcome-equation rank check now includes the constant that the cutpoints absorb: a set of columns that sums to a constant (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).
• Bootstrap relabelling of duplicated units uses a suffix guaranteed absent from the original labels; 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).

Category-specific split equations, tests, scores, bootstrap (2026-08-22)

• 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.

Inference and comparison

• 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().

Quantities of interest and diagnostics

• 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().

Data

• 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).

Documentation

• Four vignettes: 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).
• A package overview page (?"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.

Robustness (adversarial pass)

• 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.

Validation

• Bivariate-normal probabilities agree with 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.