The problem
Ordered outcomes in the social sciences often have one category that is too full: a bottom category of “no violence” that holds both structurally peaceful countries and countries at risk that happened to stay calm; a middle “neither agree nor disagree” that holds genuine neutrals and face-saving non-responders; a top category of “maximum provisions” that holds step-by-step adopters and members who adopted everything at once for other reasons. A standard ordered probit or logit treats all of these as one homogeneous category and, when the mixture is real, delivers biased coefficients and misleading predicted probabilities.
The inflated ordered models of Harris and Zhao (2007), Bagozzi and
Mukherjee (2012), and Bagozzi, Joo and Mukherjee (2024) add a second
latent equation: a split (inflation) equation that
sorts units into the ordered regime or the inflated regime, where the
inflated category is observed with certainty. iop fits
these models for any single inflated category, with probit
(iop()) or logit (iol()) links, optional
correlated errors, random intercepts, and the same methods as the plain
ordered models (oprobit(), ologit()) it also
provides.
Four estimators, one interface
| Function | Model | Extras |
|---|---|---|
oprobit() |
ordered probit |
parallel = (partial proportional odds) |
ologit() |
ordered logit | parallel = |
iop() |
inflated ordered probit (any category) |
correlated = TRUE (ZiOPC/MiOPC/TiOPC) |
iol() |
inflated ordered logit (any category) |
All four accept the same formula, y ~ x1 + x2 | z1 + z2,
where the part after the | is the inflation equation (it
always has an intercept), plus weights,
offset, re = (random intercepts),
fe = (unit fixed effects), and
se = "analytic" | "robust" | "cluster". They return one
class, "iord", so every method below works identically on
every model.
Zero inflation: political violence
The bp data (Besley and Persson 2009, as analyzed by
Bagozzi, Hill, Moore and Mukherjee 2015) record whether a country-year
saw no political violence, repression, or civil war.
library(iop)
data(bp)
round(prop.table(table(bp$violence)), 3)
#>
#> none repression civil war
#> 0.716 0.195 0.089The bottom category holds 72 percent of country-years. Fit the plain ordered probit and the zero-inflated ordered probit (the published specification uses the same covariates in both equations):
f_op <- violence ~ loggdppc + parliament + disaster + major_oil + major_primary
f_ziop <- violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
loggdppc + parliament + disaster + major_oil + major_primary
m_op <- oprobit(f_op, data = bp)
m_ziop <- iop(f_ziop, data = bp, inflate = "bottom")
summary(m_ziop)
#>
#> Inflated ordered probit (inflated category: none)
#> Call: iop(formula = f_ziop, data = bp, inflate = "bottom")
#> Response levels (in order): none < repression < civil war
#> N = 1984 inference: analytic
#>
#> Outcome equation (ordered probit):
#> Estimate Std. Error z value Pr(>|z|)
#> loggdppc 0.001604 0.057695 0.028 0.9778
#> parliament -0.095934 0.168983 -0.568 0.5702
#> disaster 0.273732 0.032817 8.341 < 2e-16 ***
#> major_oil 1.900762 0.450193 4.222 2.42e-05 ***
#> major_primary -0.562976 0.255867 -2.200 0.0278 *
#>
#> Cutpoints:
#> Estimate Std. Error z value Pr(>|z|)
#> none|repression 0.5017 0.4151 1.209 0.226854
#> repression|civil war 1.3985 0.4244 3.296 0.000982 ***
#>
#> Inflation equation (P(ordered regime); inflated category "none", observed share 0.716, mean fitted P(ordered regime) 0.793):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 20.6179 4.1970 4.913 8.99e-07 ***
#> loggdppc -2.2461 0.4425 -5.076 3.85e-07 ***
#> parliament -0.4388 0.4236 -1.036 0.300
#> disaster -0.1591 0.1700 -0.936 0.349
#> major_oil -4.7601 42.0348 -0.113 0.910
#> major_primary 3.9884 42.0309 0.095 0.924
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> logLik = -1384.26 AIC = 2794.51 BIC = 2867.22 df = 13The inflation coefficients predict the probability of the
ordered regime – positive means more likely to be a
country at risk of violence, less likely to be structurally peaceful –
and predict(type = "inflated") returns the complementary
probability. The split sharpens the outcome equation: GDP per capita, a
strong predictor in the plain model, now works through the split (richer
countries are more likely to be structurally peaceful) rather than
through the ordered stage.
Is the inflation stage needed?
The plain model is the limit of the inflated one as the inflation
intercept goes to infinity – a point outside the parameter space – so
the usual likelihood-ratio test does not apply.
inflation_test() refits the plain model internally and
reports the Vuong test with its AIC and BIC corrections, together with
the information criteria:
inflation_test(m_ziop)
#> Inflation test: inflated category "none" (observed share 0.716)
#>
#> model logLik df AIC BIC
#> Inflated ordered probit (inflated category: none) -1384.26 13 2794.51 2867.22
#> Ordered probit (uninflated) -1432.24 7 2878.48 2917.63
#>
#> Vuong test, inflated vs uninflated (positive favors the inflated model):
#> correction statistic p_inflated_better p_ordered_better p_two_sided
#> raw 5.256 0.000 1.000 0.000
#> AIC 4.599 0.000 1.000 0.000
#> BIC 2.761 0.003 0.997 0.006
#> (The plain model is nested at the boundary of the inflated one, so the Vuong test is contested for this comparison; for a reported test add the parametric-bootstrap LR with boot = 199 or more.)
compare_models(op = m_op, ziop = m_ziop)
#> model type logLik df AIC
#> 1 op Ordered probit -1432.241 7 2878.483
#> 2 ziop Inflated ordered probit (inflated category: none) -1384.256 13 2794.511
#> BIC N
#> 1 2917.633 1984
#> 2 2867.219 1984Quantities of interest by regime
Predicted probabilities are available for each category, for the regime, and for the posterior probability that an observed zero is a structural (inflated) zero:
head(predict(m_ziop, type = "prob"), 3)
#> none repression civil war
#> [1,] 0.9674842 0.022065503 0.010450287
#> [2,] 0.9822735 0.008319777 0.009406761
#> [3,] 0.9769820 0.015619602 0.007398424
summary(predict(m_ziop, type = "inflated"))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000e+00 4.400e-07 5.295e-03 2.073e-01 2.229e-01 1.000e+00
summary(predict(m_ziop, type = "posterior")[bp$violence == "none"])
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000e+00 7.300e-07 3.021e-02 2.896e-01 6.880e-01 1.000e+00First differences move one covariate in the outcome equation, the inflation equation, or both, and report the change in every category probability and in the regime probability with delta-method intervals:
first_difference(m_ziop, "loggdppc", from = 7, to = 9)
#> First differences: loggdppc from 7 to 9 [at the typical profile]
#> component from to diff lower upper method
#> P(y = none) 0.645 0.781 0.136 -0.080 0.351 delta
#> P(y = repression) 0.253 0.156 -0.097 -0.251 0.058 delta
#> P(y = civil war) 0.102 0.063 -0.039 -0.102 0.024 delta
#> P(ordered regime) 1.000 0.616 -0.384 -1.012 0.243 delta
first_difference(m_ziop, "loggdppc", from = 7, to = 9, stage = "inflation")
#> First differences: loggdppc from 7 to 9 (inflation equation only) [at the typical profile]
#> component from to diff lower upper method
#> P(y = none) 0.645 0.782 0.136 -0.088 0.360 delta
#> P(y = repression) 0.253 0.156 -0.097 -0.257 0.063 delta
#> P(y = civil war) 0.102 0.063 -0.039 -0.104 0.025 delta
#> P(ordered regime) 1.000 0.616 -0.384 -1.012 0.243 deltavignette("quantities") covers predictions on covariate
profiles, standard errors for predicted probabilities, average marginal
effects, plots, and regression tables.
The logit versions
ologit() and iol() are the logit
counterparts, with the same formula and methods (there is no
correlated-errors option for the logit link, which has no canonical
bivariate form). Because all four models share one engine, their fit
statistics are directly comparable:
m_ol <- ologit(f_op, data = bp)
m_ziol <- iol(f_ziop, data = bp, inflate = "bottom")
compare_models(op = m_op, ziop = m_ziop, ol = m_ol, ziol = m_ziol)
#> model type logLik df AIC
#> 1 op Ordered probit -1432.241 7 2878.483
#> 2 ziop Inflated ordered probit (inflated category: none) -1384.256 13 2794.511
#> 3 ol Ordered logit -1435.443 7 2884.885
#> 4 ziol Inflated ordered logit (inflated category: none) -1383.536 13 2793.071
#> BIC N
#> 1 2917.633 1984
#> 2 2867.219 1984
#> 3 2924.035 1984
#> 4 2865.778 1984Top inflation: escape-flexibility provisions
Any category can be inflated. The pta data (Baccini, Dur
and Elsig 2015, as analyzed by Bagozzi, Joo and Mukherjee 2024) count
the escape-flexibility provisions in 559 preferential trade agreements;
the top category – all four provisions – holds 46 percent of agreements
and mixes step-by-step insurance-seekers with members that adopt maximum
flexibility at once:
data(pta)
table(pta$flexibility)
#>
#> 0 1 2 3 4
#> 66 66 80 92 255
m_tiop <- iop(flexibility ~ depth * democracy + gdp + gdppc + trade + gattwto + members +
democratization | gdp + gdppc + democracy + democratization,
data = pta, inflate = "top")
summary(m_tiop)
#>
#> Inflated ordered probit (inflated category: 4)
#> Call: iop(formula = flexibility ~ depth * democracy + gdp + gdppc +
#> trade + gattwto + members + democratization | gdp + gdppc +
#> democracy + democratization, data = pta, inflate = "top")
#> Response levels (in order): 0 < 1 < 2 < 3 < 4
#> N = 559 inference: analytic
#>
#> Outcome equation (ordered probit):
#> Estimate Std. Error z value Pr(>|z|)
#> depth 0.076095 0.011150 6.824 8.83e-12 ***
#> democracy 0.233859 0.202387 1.156 0.24788
#> gdp 0.123126 0.054717 2.250 0.02443 *
#> gdppc -0.198259 0.093155 -2.128 0.03332 *
#> trade -0.050652 0.029668 -1.707 0.08776 .
#> gattwto 0.222431 0.141240 1.575 0.11529
#> members -0.001499 0.007050 -0.213 0.83165
#> democratization -0.386944 0.161517 -2.396 0.01659 *
#> depth:democracy -0.032583 0.012046 -2.705 0.00683 **
#>
#> Cutpoints:
#> Estimate Std. Error z value Pr(>|z|)
#> 0|1 0.1440 1.0281 0.140 0.8886
#> 1|2 0.7472 1.0275 0.727 0.4671
#> 2|3 1.3770 1.0282 1.339 0.1805
#> 3|4 2.1091 1.0363 2.035 0.0418 *
#>
#> Inflation equation (P(ordered regime); inflated category "4", observed share 0.456, mean fitted P(ordered regime) 0.764):
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.9215 1.3383 -0.689 0.491090
#> gdp 0.4359 0.1192 3.658 0.000254 ***
#> gdppc -0.7506 0.2282 -3.290 0.001002 **
#> democracy -1.3982 0.3824 -3.657 0.000256 ***
#> democratization -0.7269 0.3531 -2.059 0.039512 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> logLik = -704.11 AIC = 1444.21 BIC = 1522.08 df = 18Here the inflation equation contains covariates that are also in the
outcome equation, but the outcome equation has covariates the inflation
equation lacks, so the split is not identified by functional form alone.
The repression data (Girod, Stewart and Walters 2018, the
appendix application of the same article) are a second top-inflated
example:
data(repression)
m_rep <- iop(repression ~ negxpol * oilrent | dom_media + civil_war,
data = repression, inflate = "top")
inflation_test(m_rep)$vuong
#> correction statistic p_inflated_better p_ordered_better p_two_sided
#> 1 raw 3.625263 0.0001443336 0.9998557 0.0002886672
#> 2 AIC 3.286201 0.0005077434 0.9994923 0.0010154868
#> 3 BIC 2.624119 0.0043436687 0.9956563 0.0086873375Middle inflation: a survey example
The middle-inflated model of Bagozzi and Mukherjee (2012) treats the
middle category of an attitude scale as a mixture of genuine neutrals
and non-responders. The Eurobarometer data of that article cannot be
redistributed, so here is a simulated version with riop(),
which draws from the package’s own data-generating process:
set.seed(1)
d <- riop(1000, beta = c(0.8, -0.5), tau = c(-0.8, 0.8), gamma = c(0.2, 1),
inflate = "middle", labels = c("disagree", "neutral", "agree"))
table(d$y)
#>
#> disagree neutral agree
#> 148 706 146
m_miop <- iop(y ~ x1 + x2 | z1, data = d, inflate = "middle")
m_miop
#> Inflated ordered probit (inflated category: neutral)
#> Response levels (in order): disagree < neutral < agree
#> Call: iop(formula = y ~ x1 + x2 | z1, data = d, inflate = "middle")
#>
#> Outcome coefficients:
#> x1 x2
#> 0.8521 -0.6618
#>
#> Cutpoints:
#> disagree|neutral neutral|agree
#> -0.8214 0.8467
#>
#> Inflation coefficients (positive = higher P(ordered regime)):
#> (Intercept) z1
#> 0.0285 0.8739
#>
#> logLik: -605.13 N: 1000The true values are beta = (0.8, -0.5), cutpoints
(-0.8, 0.8), and inflation coefficients
(0.2, 1). The posterior probability that an observed
“neutral” is a non-response rather than a genuine neutral is
predict(type = "posterior"):
summary(predict(m_miop, type = "posterior")[d$y == "neutral"])
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.0299 0.5208 0.7534 0.6898 0.9033 0.9994inflate takes "bottom",
"middle" (odd number of categories), "top", or
a category label (inflate = "neutral" here), so any single
category can be inflated.
Beyond political science: two examples from economics and health
Nothing in the models is specific to political science; the bundled data sets are simply the ones behind the articles that introduced them. Two standard data sets from the AER package show the same workflow on an economics and a health outcome (the chunks run only if AER is installed).
The Fair (1978) extramarital-affairs data (AER::Affairs,
n = 601) record the number of affairs in the past year – 0, 1, 2, 3, 7,
or 12, an ordinal response with three quarters of the sample at zero. A
zero-inflated ordered probit separates the decision to have any affair
(the inflation equation, here driven by gender and children, which are
excluded from the outcome equation) from how many, given the propensity
(the ordered equation):
data("Affairs", package = "AER")
m_aff <- iop(affairs ~ age + yearsmarried + religiousness + rating | gender + children,
data = Affairs, inflate = "bottom")
m_aff
#> Inflated ordered probit (inflated category: 0)
#> Response levels (in order): 0 < 1 < 2 < 3 < 7 < 12
#> Call: iop(formula = affairs ~ age + yearsmarried + religiousness +
#> rating | gender + children, data = Affairs, inflate = "bottom")
#>
#> Outcome coefficients:
#> age yearsmarried religiousness rating
#> -0.0305 0.0895 -0.2683 -0.3478
#>
#> Cutpoints:
#> 0|1 1|2 2|3 3|7 7|12
#> -2.1370 -1.8465 -1.6841 -1.4892 -0.9205
#>
#> Inflation coefficients (positive = higher P(ordered regime)):
#> (Intercept) gendermale childrenyes
#> 0.0024 0.4422 0.2282
#>
#> logLik: -527.84 N: 601
colMeans(predict(m_aff, type = "zeros")) # zeros from the inflation process / the ordered stage
#> inflation ordered
#> 0.3576842 0.3910722
inflation_test(m_aff)
#> Inflation test: inflated category "0" (observed share 0.750)
#>
#> model logLik df AIC BIC
#> Inflated ordered probit (inflated category: 0) -527.84 12 1079.69 1132.47
#> Ordered probit (uninflated) -531.21 9 1080.43 1120.02
#>
#> Vuong test, inflated vs uninflated (positive favors the inflated model):
#> correction statistic p_inflated_better p_ordered_better p_two_sided
#> raw 1.210 0.113 0.887 0.226
#> AIC 0.133 0.447 0.553 0.894
#> BIC -2.236 0.987 0.013 0.025
#> (The plain model is nested at the boundary of the inflated one, so the Vuong test is contested for this comparison; for a reported test add the parametric-bootstrap LR with boot = 199 or more.)The two types of zeros are roughly equally common here, and the
inflated model improves the fit only modestly over the plain ordered
probit: the Vuong statistic leans toward it, the AIC is close, and the
BIC prefers the simpler model – a reminder that
inflation_test(boot = 199) (the bootstrap likelihood ratio,
not run here) is the test to report.
Self-rated health in the 1988 National Medical Expenditure Survey
(AER::NMES1988, n = 4,406) is poor / average / excellent,
with four fifths of respondents at “average” – the middle category of a
survey scale where “average” can also be a default answer. The
middle-inflated model lets education and age move the probability of a
considered answer, and the health covariates move the answer itself:
data("NMES1988", package = "AER")
m_nmes <- iop(health ~ age + gender + married + chronic + adl + income + insurance | school + age,
data = NMES1988, inflate = "middle")
round(summary(m_nmes)$coefficients[, 1:2], 3)
#> Estimate Std. Error
#> age -0.018 0.042
#> gendermale 0.012 0.053
#> marriedyes -0.152 0.057
#> chronic -0.392 0.030
#> adllimited -0.932 0.080
#> income 0.050 0.012
#> insuranceyes 0.273 0.061
#> poor|average -1.789 0.320
#> average|excellent 0.781 0.326
#> infl_(Intercept) 0.923 0.721
#> infl_school -0.032 0.013
#> infl_age -0.021 0.082
ame(m_nmes, vars = c("chronic", "insurance"))
#> Average marginal effects (delta-method intervals; derivatives for continuous covariates, discrete changes otherwise)
#> variable contrast component estimate lower upper
#> chronic dP/dx P(y = poor) 0.0518 0.0464 0.0571
#> chronic dP/dx P(y = average) -0.0065 -0.0112 -0.0019
#> chronic dP/dx P(y = excellent) -0.0452 -0.0508 -0.0397
#> chronic dP/dx P(ordered regime) 0.0000 0.0000 0.0000
#> insurance no -> yes P(y = poor) -0.0383 -0.0558 -0.0208
#> insurance no -> yes P(y = average) 0.0093 0.0024 0.0161
#> insurance no -> yes P(y = excellent) 0.0290 0.0175 0.0406
#> insurance no -> yes P(ordered regime) 0.0000 0.0000 0.0000Everything downstream – predict(),
first_difference(), ame(),
classification(), the tests of the previous sections –
applies unchanged.
Where next
-
vignette("quantities")– predicted probabilities with standard errors, the two-types-of-zeros decomposition (predict(type = "zeros"),decompose = TRUE), first differences, average marginal effects, plots, regression tables (broom, modelsummary, texreg), classification and accuracy scores (classification()), and simulated-residual diagnostics. -
vignette("panels")– random intercepts, unit fixed effects and the split-panel jackknife, the Mundlak device, cluster-robust and bootstrap standard errors, and the package’s Monte Carlo on short panels. -
vignette("model")– the likelihood, the sign conventions, the category-specific split equations ofsplit = "category"andsplit_test(), identification and exclusion restrictions, multi-start estimation and boundary cases, which test for which comparison (Vuong, bootstrap likelihood ratio, LM/LR), partial proportional odds, the names these models go by elsewhere, and how the package was validated.
References
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.
Baccini, L., Dur, A. and Elsig, M. (2015). The politics of trade agreement design: Revisiting the depth-flexibility nexus. International Studies Quarterly, 59, 765-775.
Besley, T. and Persson, T. (2009). Repression or civil war? American Economic Review: Papers and Proceedings, 99, 292-297.
Girod, D.M., Stewart, M.A. and Walters, M.R. (2018). Mass protests and the resource curse: The politics of demobilization in rentier autocracies. Conflict Management and Peace Science, 35, 503-522.
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.