This vignette is the reference for what iop estimates
and how: the likelihood and its sign conventions, what identifies the
split, how the optimizer guards against the pitfalls of mixture
likelihoods, what the boundary warnings mean, which test to use for
which comparison, and how the implementation was validated.
For an outcome with categories \(j = 0, \ldots, J-1\) and covariates \(x_i\), the latent equation is \(y_i^* = x_i'\beta + \epsilon_i\), with \(y_i = j\) when \(\tau_{j-1} < y_i^* \le \tau_j\) (\(\tau_{-1} = -\infty\), \(\tau_{J-1} = \infty\)). With \(F\) the standard normal (probit) or logistic (logit) distribution function, the cell probabilities are
\[ \pi_{ij} = F(\tau_j - x_i'\beta) - F(\tau_{j-1} - x_i'\beta). \]
The outcome equation has no intercept – the cutpoints absorb it –
matching MASS::polr() and ordinal::clm().
Internally the cutpoints are parameterized as \(\tau_0\) and the logs of their increments
whenever every term is parallel, so they stay ordered during
optimization; the reported cutpoints and their covariance are on the
natural scale.
parallel = FALSE, or a one-sided formula naming the
terms held parallel (parallel = ~ . - x2 relaxes only
x2), gives the relaxed terms a separate coefficient per
cutpoint, \(P(y_i \le j) = F(\tau_j -
x_i'\beta - \tilde x_i'\beta_j)\): the partial
proportional-odds / generalized ordered model of Peterson and Harrell
(1990) and Williams (2006). parallel_test() reports the
likelihood-ratio test of each term’s parallel restriction (the
likelihood counterpart of the Brant test) and the omnibus test, and
parallel = "auto" relaxes terms one at a time while the
smallest p-value is below 0.05 (a forward version of Stata’s
gologit2, autofit), storing the steps in
$autofit:
library(iop)
data(pta)
m_pl <- ologit(flexibility ~ depth + democracy + gdp + gdppc + trade, data = pta)
parallel_test(m_pl)
#> Likelihood-ratio tests of the parallel-regression assumption (relax one term at a time; last row relaxes all)
#> term LR df p.value logLik
#> depth 54.021 3 1.11e-11 -702.60
#> democracy 8.041 3 4.52e-02 -725.59
#> gdp 19.753 3 1.91e-04 -719.73
#> gdppc 9.552 3 2.28e-02 -724.83
#> trade 12.223 3 6.66e-03 -723.49
#> all terms 80.234 15 6.33e-11 -689.49
m_auto <- ologit(flexibility ~ depth + democracy + gdp + gdppc + trade, data = pta, parallel = "auto")
m_auto$autofit$relaxed
#> [1] "depth" "democracy"
lr_test(m_pl, m_auto)
#> Likelihood-ratio test
#> LR = 66.0039 on 6 df, p = 2.69e-12Non-parallel fits can imply negative cell probabilities where the category-specific curves cross; the fit warns when that happens at the optimum, and the remedy is to hold more terms parallel.
A second latent equation, \(s_i^* = z_i'\gamma + u_i\), assigns unit \(i\) to the ordered regime (\(s_i = 1\) when \(s_i^* > 0\), with probability \(P(s_i = 1) = F(z_i'\gamma)\)) or to the inflated regime (\(s_i = 0\)), in which the outcome is the inflated category \(k\) with certainty. Hence
\[ P(y_i = j) = F(z_i'\gamma)\,\pi_{ij} + 1\{j = k\}\,[1 - F(z_i'\gamma)], \]
the zero-inflated ordered probit of Harris and Zhao (2007) for \(k = 0\), the middle-inflated model of
Bagozzi and Mukherjee (2012), the top-inflated model of Bagozzi, Joo and
Mukherjee (2024), and in general the inflated ordered outcome of Brooks,
Harris and Spencer (2012). Sign convention: following
Harris and Zhao (2007) and the political-science literature, positive
inflation coefficients raise the probability of the ordered
(non-inflated) regime; predict(type = "inflated") returns
\(1 - F(z_i'\gamma)\).
The inflated category is observed from two sources, so the posterior
probability that an observation in category \(k\) is an inflated case is \(P(s_i = 0 \mid y_i = k) = [1 - F(z_i'\gamma)]
/ P(y_i = k)\) (predict(type = "posterior")).
The models above use one split equation for every unit. Brown, Harris and Spencer (2020) generalise this: every non-inflated category \(j\) gets its own split equation \(s_{ij}^* = z_i'\gamma_j + u_{ij}\), so a unit whose ordered outcome would be \(j\) is “tempered” into the inflated category with its own probability \(1 - F(z_i'\gamma_j)\),
\[ P(y_i = j) = F(z_i'\gamma_j)\,\pi_{ij} \ (j \ne k), \qquad P(y_i = k) = 1 - \sum_{j \ne k} F(z_i'\gamma_j)\,\pi_{ij}, \]
with a \(\rho_j\) per equation in
the correlated version. This is the generalised zero-/middle-inflated
ordered probit (GZiOP, GMiOP, and the correlated GZiOPC/GMiOPC); the
common-split model is the special case in which all \(\gamma_j\) (and \(\rho_j\)) are equal.
split = "category" fits it, starting from the common-split
fit (kept in $loglik_common) and checking that the
likelihood is not below it; split_test() tests the
restriction both by a Lagrange-multiplier (score) test computed from the
common-split fit alone and by the likelihood-ratio test, which is
standard here because the restriction is interior. Brown, Harris and
Spencer find the generalised model favoured over the original ZIOP on
Harris and Zhao’s tobacco data and over the MIOP on the Eurobarometer
data of Bagozzi and Mukherjee (2012). The Besley–Persson
political-violence data behave the same way. With the published
two-covariate split (log GDP per capita and parliamentary
democracy):
f_bp <- violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
loggdppc + parliament
m_common <- iop(f_bp, data = bp, inflate = "bottom")
m_category <- iop(f_bp, data = bp, inflate = "bottom", split = "category")
#> Warning in .iord_fit(formula, data, link = "probit", inflate = inflate, : iop:
#> the split equation for category civil war is quasi-separated (standardized
#> coefficients implausibly large; the split acts as a step function and its
#> estimate may not be finite) -- consider split = "common".
compare_models(common = m_common, category = m_category)
#> model
#> 1 common
#> 2 category
#> type
#> 1 Inflated ordered probit (inflated category: none)
#> 2 Inflated ordered probit (inflated category: none), category-specific split
#> logLik df AIC BIC N
#> 1 -1385.909 10 2791.818 2847.747 1984
#> 2 -1357.391 13 2740.782 2813.489 1984
split_test(m_common)
#> Test of a common split equation against category-specific split equations (Brown, Harris and Spencer 2020)
#> LM information: outer product of gradients (Brown, Harris and Spencer); a local screening statistic -- report the LR when they disagree
#> test statistic df p.value logLik_common logLik_category
#> LM 266.797 3 1.52e-57 -1385.91 NA
#> LR 57.036 3 2.52e-12 -1385.91 -1357.39
m_category
#> Inflated ordered probit (inflated category: none), category-specific split
#> Response levels (in order): none < repression < civil war
#> Call: iop(formula = f_bp, data = bp, inflate = "bottom", split = "category")
#>
#> Outcome coefficients:
#> loggdppc parliament disaster major_oil major_primary
#> 0.1342 -0.0647 0.2582 1.9012 -0.4550
#>
#> Cutpoints:
#> none|repression repression|civil war
#> 1.4492 2.2816
#>
#> Inflation coefficients (positive = higher P(ordered regime); one split equation per non-inflated category):
#> repression civil war
#> (Intercept) 16.1068 103.3873
#> loggdppc -1.7225 -12.1321
#> parliament -0.9057 3.1979
#>
#> logLik: -1357.39 N: 1984
#> WARNING: boundary fit -- the split equation for category civil war is quasi-separated (standardized coefficients implausibly large; the split acts as a step function and its estimate may not be finite) -- consider split = "common".Both tests reject the common split decisively. The LM is a
local test evaluated at the common-split estimates – correctly
sized under the null in simulation (with either the outer-product
information of Brown, Harris and Spencer, the default, or the observed
information, information = "hessian"), but far from the LR
when the restriction is strongly violated, as here (267 versus 57; the
observed information is not even positive definite at that point, so the
Hessian form is reported as NA). Treat the LM as a screen
and report the LR. The fitted split equations say why the common split
is rejected: the split that governs whether a country can experience
repression is the familiar gradual one, while the split that
governs whether it can experience civil war is far steeper in
income – civil war is effectively ruled out above a threshold of log GDP
per capita, an empirical regularity the common split had to average
over. Under a category-specific split the coefficients are named
infl_<term>:<category>, and
predict(type = "regime") and the regime rows of
first_difference() / ame() report one
probability per non-inflated category.
Two cautions, both from Brown, Harris and Spencer. Identification now
rests on \(J - 1\) split equations, so
exclusion restrictions and the number of observations in the rarer
categories matter more: with the full five-covariate split on these data
the civil-war equation, fed by 176 civil-war years, acquires a
near-cancelling pair of coefficients on the two exporter dummies and the
optimizer reports an unclean optimum – the signature of a split the data
cannot pin down, which the convergence and boundary warnings surface.
And a large gap between the LM and LR statistics is itself a warning
about the category-specific fit. Random intercepts are not available
with split = "category".
The split is identified by the shape of the two link functions and,
far more credibly, by an exclusion restriction: a
covariate in one equation but not the other. Two advisory messages guard
this. If the formula has no | part, the inflation equation
reuses the outcome covariates; if the | part contains no
covariate absent from the outcome equation, the split rests on
functional form alone. Both cases fit, but say so:
data(bp)
m_same <- iop(violence ~ loggdppc + disaster | loggdppc + disaster, data = bp, inflate = "bottom")
#> The inflation equation contains no covariate that is excluded from the outcome equation; the split is then identified by functional form alone. An exclusion restriction is advisable.The published zero-inflated ordered probit of political violence uses
the same covariates in both equations, so vignette("iop")
triggers this message deliberately; when the design offers a genuine
exclusion, use it.
The log-likelihood, its analytic gradient, and the per-observation scores are implemented in C++ for every combination of link, inflation, correlation, parallel / non-parallel terms, and random intercepts. Every fit column-scales the design, runs BFGS with the analytic gradient, and finishes with exact-Hessian Newton steps (numerical Jacobian of the analytic gradient) with step halving; convergence is declared from the Newton decrement \(g'H^{-1}g\), so the optimum is located to machine precision even when quasi-Newton progress along a flat ridge has stalled.
The likelihood of an inflated mixture can have several local maxima:
a “soft” split with modest inflation coefficients and a “sharp” split
with steep ones can both be stationary points. Every inflated fit is
therefore a multi-start. The inflation equation is started from a binary
model of membership in the non-inflated category at several slope scales
and signs, from flat high-regime values, and from the plain ordered
baseline; each start gets a short quasi-Newton run, the two best are run
to convergence and polished, and the best optimum is kept.
$start_logliks records the log-likelihood reached from
every start (short-run values for the non-finalists). The Besley–Persson
zero-inflated ordered probit with the full published specification is a
real example of two modes:
m_ziop <- iop(violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
loggdppc + parliament + disaster + major_oil + major_primary,
data = bp, inflate = "bottom")
sort(round(m_ziop$start_logliks, 2), decreasing = TRUE)
#> [1] -1384.26 -1384.26 -1390.35 -1390.35 -1390.35 -1390.38 -1390.39 -1401.75
#> [9] -1412.71Several starts reach a local optimum near -1390; the reported fit is
the mode at -1384.26 (Table 1 of Bagozzi et al. 2015). A spread of
several units across starts is the signature of multimodality;
start = lets you add starts of your own.
Every inflated fit also fits its plain ordered counterpart first
(whose log-likelihood is concave) and checks that the inflated fit is
not below it: the plain model is the limit of the inflated one as the
inflation intercept tends to infinity, so the inflated log-likelihood
can never be lower at a true optimum. The baseline is kept in
$loglik_uninflated and is what
inflation_test() compares against.
Some data contain no identifiable inflation process. The fit then drifts toward the no-inflation limit – every regime probability at 1, the inflated log-likelihood equal to the ordered baseline – and is flagged:
set.seed(1)
d0 <- riop(500, beta = c(0.8, -0.5), tau = c(-0.5, 0.7), gamma = c(3.5, 0.1), inflate = "bottom")
mean(attr(d0, "regime") == 0) # no inflated-regime units were drawn
#> [1] 0
m0 <- iop(y ~ x1 + x2 | z1, data = d0, inflate = "bottom")
#> Warning in .iord_fit(formula, data, link = "probit", inflate = inflate, : iop:
#> the inflation equation is degenerate (regime probabilities at 1 for every
#> unit): the data show no separate inflation process and the fit collapses toward
#> the plain ordered model.
m0$boundary
#> [1] TRUE
c(inflated = m0$loglik, ordered = m0$loglik_uninflated)
#> inflated ordered
#> -438.8916 -438.8916A second boundary has the opposite look: the split acts as a step
function – a split covariate that classifies units almost perfectly, or
a few observations in the inflated category predicted to be inflated
with certainty – so the split coefficients run off to infinity
(standardized coefficients in the tens or hundreds, standard errors in
the hundreds or, once the curvature collapses, numerically zero). The
flag catches this graded case too, by the size of the standardized split
coefficients (slope above 10, intercept above 40; the sharpest
legitimate splits in the package’s applications reach 6 and 21) and by
degenerate split standard errors (above 50 on the standardized scale),
and the remedy is a simpler inflation equation. The same signatures
(standardized coefficient above 10, standardized standard error above
50) are applied to the outcome equation, where a covariate that
classifies an extreme category almost perfectly has the same non-finite
estimate (unit dummies are exempt). The same $boundary flag
marks a correlation estimate at \(|\rho| >
0.985\) – report the uncorrelated fit – and a random-intercept
standard deviation at zero – there is no unit heterogeneity in that
equation.
A third boundary belongs to the cutpoints rather than the split. The
ordered stage can empty the inflated category: the inflation process
absorbs every observation in it, the ordered-stage probability of that
category goes to zero for every unit, and the fit is a hurdle-type split
– the adjacent cutpoint runs to \(-\infty\) (bottom-inflated) or \(+\infty\) (top-inflated), or the two
cutpoints bracketing a middle category coincide. The likelihood is flat
in that direction, the finite cutpoint the optimizer stops at is an
artifact of the stopping rule, and a generalized inverse would hand it a
reassuring standard error. The flag fires when the ordered stage
supplies less than 1e-6 of the inflated category’s fitted probability
mass (predict(type = "zeros") shows the two components), or
less than 1e-2 of it while the adjacent cutpoint’s standard error
already exceeds 50 (the optimizer stopped on its gradient tolerance
before the cutpoint ran further out); it names the cutpoint and reports
its standard error as NA; the other parameters are
unaffected. Our stress battery found this case on data with no genuine
inflation (a top-inflated fit to AER::BankWages, where
every manager is attributed to the split), which is the situation in
which it should be expected.
Two further fields complete the picture and are printed by
summary(): $ill_conditioned (the information
matrix has reciprocal condition number below 1e-12 in the optimizer’s
scaled parameterization, with every covariate column standardized so
that the units of the covariates play no role – the Hessian is
differenced in that parameterization too, so the standard errors are
likewise invariant to a change of units: near-collinear covariates or a
near-flat split direction; those standard errors are unreliable) and
$flat_hessian (a converged optimum with a near-flat
direction: weak identification, typically an inflation equation without
an exclusion restriction). A delta-method variance that is not positive
at such an optimum is reported as an NA standard error,
never as zero. All of these thresholds are fixed by design and listed in
?iop (“Diagnostics and fixed thresholds”); they were set so
that every legitimate fit in the package’s applications and stress tests
passes and every constructed boundary case is caught.
| Comparison | Tool | Why |
|---|---|---|
| inflated vs plain ordered | inflation_test() (Vuong, AIC/BIC, bootstrap LR) |
the restriction sits on the boundary; no standard LR distribution |
| correlated vs uncorrelated inflated probit | lr_test() |
\(\rho = 0\) is interior; one degree of freedom |
| common vs category-specific split | split_test() (LM and LR) |
interior restriction (equal split coefficients) |
| non-parallel vs parallel terms | parallel_test(), lr_test() |
interior restriction |
| different inflated categories, probit vs logit | vuong(), compare_models(),
classification() |
non-nested |
vuong() reports the raw Vuong (1989) statistic and the
AIC- and BIC-corrected versions, as pscl::vuong() does;
lr_test() refuses pairs that differ in link, inflated
category, or random-intercept structure.
The inflated-versus-plain comparison deserves a word. The applied
literature (Harris and Zhao 2007; Bagozzi et al. 2015) uses the Vuong
test, and so does inflation_test(). Wilson (2015) and Dale
and Sirchenko (2021) object that the plain model is nested in
the inflated one (at the boundary), while the Vuong z-test is derived
for non-nested or overlapping models; Stata’s ziop2 offers
no Vuong option for this pair, and Dale and Sirchenko’s Monte Carlo
finds the information criteria the most reliable selectors. The direct
remedy for a boundary restriction is a parametric bootstrap of the
likelihood-ratio statistic (Andrews 2001): simulate from the fitted
plain ordered model, refit both models, and take the share of bootstrap
statistics at least as large as the observed one.
inflation_test(boot = 199) does this; it costs 199 inflated
refits (a few minutes for bp), so it is opt-in:
Report all three – Vuong, information criteria, bootstrap LR – when
the case for inflation is the point of the analysis.
classification() adds the proper scoring rules (Brier,
ranked probability) and the classification table for any pair of fits,
in or out of sample.
weights enter the log-likelihood as frequency weights: a
fit with weight 2 equals a fit on duplicated rows, in the estimates, the
model-based standard errors, the information criteria (whose \(n\) is the weight total), the Vuong
statistics, and the averaged quantities of interest. Survey (sampling,
probability) weights are a different thing: with them the default
model-based standard errors are not the design-based ones, so use
se = "robust" (the pseudo-maximum-likelihood sandwich) or
cluster = for the primary sampling units, as one would
after pweight in Stata. offset and
offset_inflation enter the two latent equations, and
se selects analytic (inverse observed information), robust
(sandwich over observations, or over units for random-intercept fits),
cluster-robust (cluster = implies
se = "cluster"), or nonparametric bootstrap standard errors
(se = "bootstrap", nboot refits over the rows,
the clusters, or the random-intercept units, with percentile intervals
from confint(type = "percentile"); Dale and Sirchenko 2021
find the bootstrap better calibrated than the asymptotic standard errors
for the error correlation in small samples). The covariance is formed on
the internal scale and mapped to the natural scale by the delta method;
vcov(scale = "internal") returns the former, which the
delta-method quantities of interest use.
The same models travel under several names, and one similarly named model is a different thing altogether:
In iop |
Elsewhere |
|---|---|
iop(inflate = "bottom") |
ZIOP / zero-inflated ordered probit (Harris and Zhao 2007); Stata
zioprobit; ziop2 with exogenous switching
(Dale and Sirchenko 2021); “two-part” ZIOP |
iop(..., correlated = TRUE) |
ZIOPC, “ZIOP with correlated errors”, endogenous
switching (ziop2, endoswitch) |
iop(inflate = "middle") / "top" /
label |
MIOP (Bagozzi and Mukherjee 2012; Greene, Harris and Hollingsworth
2015), TIOP (Bagozzi, Joo and Mukherjee 2024), “inflated ordered
outcomes” (Brooks, Harris and Spencer 2012),
ziop2, infcat() |
iol() |
ZIOL, zero-inflated ordered logit; Stata ziologit |
iop(..., split = "category") |
GZIOP / GMIOP, “multiple inflation processes”, “tempering” (Brown, Harris and Spencer 2020) |
parallel = FALSE / ~ . - x |
generalized ordered logit/probit, partial proportional odds
(gologit2); distinct from GZIOP, which neither nests nor is
nested by it |
| not the same model | zero-inflated bivariate ordered probit (Gurmu and Dagne 2012): two ordered outcomes modelled jointly, with a common point mass at (0, 0) – “bivariate” refers to two responses, not to the bivariate-normal errors of the ZIOPC |
The test suite checks the bivariate-normal routine against
pbivnorm and mvtnorm, every analytic gradient
against numerical differentiation,
oprobit()/ologit() against
MASS::polr, ordinal::clm, and
VGAM::vglm (including robust and cluster-robust standard
errors via sandwich), the random-intercept models against
ordinal::clmm, parameter recovery on simulated inflated
data, and – on the bundled data – the published results: Table 1 of
Bagozzi, Hill, Moore and Mukherjee (2015) on bp (ordered
probit, ZiOP, ZiOPC, and the second-mode ZiOPC2 log-likelihoods and
error correlations) and the TiOP/TiOPC log-likelihoods of Bagozzi, Joo
and Mukherjee (2024) on pta and repression.
Beyond the suite, iop() reproduces Table 1 and the
marginal-effects figures of Bagozzi and Mukherjee (2012) on the
(non-redistributable) Eurobarometer data, and the original likelihood
code of each article evaluated at iop’s estimates returns
iop’s log-likelihood to machine precision; the scripts are
in data-raw/oracles/ in the source repository.
Bagozzi, B.E. and Mukherjee, B. (2012). A mixture model for middle category inflation in ordered survey responses. Political Analysis, 20, 369-386.
Bagozzi, B.E., Hill, D.W., Moore, W.H. and Mukherjee, B. (2015). Modeling two types of peace: The zero-inflated ordered probit (ZiOP) model in conflict research. Journal of Conflict Resolution, 59, 728-752.
Bagozzi, B.E., Joo, M.M. and Mukherjee, B. (2024). Top-category inflation in ordered international relations outcomes. Foreign Policy Analysis, 20, orae006.
Andrews, D.W.K. (2001). Testing when a parameter is on the boundary of the maintained hypothesis. Econometrica, 69, 683-734.
Brooks, R., Harris, M.N. and Spencer, C. (2012). Inflated ordered outcomes. Economics Letters, 117, 683-686.
Brown, S., Harris, M.N. and Spencer, C. (2020). Modelling category inflation with multiple inflation processes: Estimation, specification, and testing. Oxford Bulletin of Economics and Statistics, 82, 1342-1361.
Dale, D. and Sirchenko, A. (2021). Estimation of nested and zero-inflated ordered probit models. Stata Journal, 21, 3-38.
Drezner, Z. and Wesolowsky, G.O. (1990). On the computation of the bivariate normal integral. Journal of Statistical Computation and Simulation, 35, 101-107.
Greene, W.H., Harris, M.N. and Hollingsworth, B. (2015). Inflated responses in measures of self-assessed health. American Journal of Health Economics, 1, 461-493.
Gurmu, S. and Dagne, G.A. (2012). Bayesian approach to zero-inflated bivariate ordered probit regression model, with an application to tobacco use. Journal of Probability and Statistics, 2012, 617678.
Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing, 14, 251-260.
Harris, M.N. and Zhao, X. (2007). A zero-inflated ordered probit model, with an application to modelling tobacco consumption. Journal of Econometrics, 141, 1073-1099.
Peterson, B. and Harrell, F.E. (1990). Partial proportional odds models for ordinal response variables. Applied Statistics, 39, 205-217.
Vuong, Q.H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica, 57, 307-333.
Williams, R. (2006). Generalized ordered logit/partial proportional odds models for ordinal dependent variables. Stata Journal, 6, 58-82.
Wilson, P. (2015). The misuse of the Vuong test for non-nested models to test for zero-inflation. Economics Letters, 127, 51-53.