Documentation site: https://bagozzib.github.io/iop/.
iop estimates ordered probit and ordered logit models
whose outcome contains an inflated category – a single
ordered category (bottom, middle, top, or any other) that mixes
observations generated by the ordered process with observations
generated by a distinct split-population process that places them in
that category regardless of the ordered mechanism. It fits the
zero-inflated ordered probit of Harris and Zhao (2007) and its middle-
and top-inflated extensions (Bagozzi and Mukherjee 2012; Bagozzi, Hill,
Moore and Mukherjee 2015; Bagozzi, Joo and Mukherjee 2024), generalized
to any inflated category and to the logit link, alongside the plain
ordered models on the same engine. It provides:
oprobit(), ologit() – the standard ordered
models, with survey weights, offsets, partial proportional-odds
(non-parallel) effects (parallel =,
parallel_test(), parallel = "auto"), and
analytic, robust, or cluster-robust standard errors;iop() – the inflated ordered probit for any single
inflated category (inflate = "bottom" | "middle" | "top" or
a category label), with optional correlated errors
(correlated = TRUE: ZiOPC / MiOPC / TiOPC) and optional
category-specific split equations
(split = "category": the generalised GZiOP / GMiOP of
Brown, Harris and Spencer 2020, tested against the common split by
split_test());iol() – the inflated ordered logit counterpart;re =, adaptive Gauss–Hermite quadrature, optionally in
both equations), unit fixed effects (fe =, with the
split-panel jackknife fe_correction = "jackknife"), the
mundlak() correlated-random-effects device, and
cluster-robust standard errors;summary(),
predict() by regime and by type of zero (with delta-method
se.fit), first_difference() and
ame() (average marginal effects, by stage – outcome
equation, inflation equation, or both – and with the two-types-of-zeros
decomposition) and their plot() methods, the
diord()/piord()/qiord()/riord()
distribution family of the response, vuong(),
lr_test(), inflation_test() (Vuong and a
parametric-bootstrap likelihood-ratio test), split_test(),
classification() (confusion table, Brier and ranked
probability scores), compare_models(),
simulate() (DHARMa-ready), ranef(), analytic /
robust / cluster-robust / bootstrap standard errors, and broom /
modelsummary / texreg integration.The likelihood, its analytic gradient, and the bivariate-normal probabilities are implemented in C++.
# development version
# install.packages("remotes")
remotes::install_github("bagozzib/iop", build_vignettes = TRUE)vignette("iop") – getting started: the four estimators
on the bundled political-violence (zero-inflated), trade-agreement
(top-inflated), and simulated survey (middle-inflated) data.vignette("quantities") – predicted probabilities with
standard errors, first differences, average marginal effects, plots,
regression tables, and simulated-residual diagnostics.vignette("panels") – random intercepts, fixed effects
and the jackknife, the Mundlak device, cluster-robust standard errors,
and the package’s Monte Carlo on short panels.vignette("model") – the likelihood and sign
conventions, identification, multi-start estimation and boundary cases,
which test for which comparison, and validation.?iop for the estimators’ full argument list;
?"iop-package" for an index.library(iop)
## political violence (none / repression / civil war): zero-inflated ordered probit
data(bp)
m_op <- oprobit(violence ~ loggdppc + parliament + disaster + major_oil + major_primary, data = bp)
m_ziop <- iop(violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
loggdppc + parliament + disaster + major_oil + major_primary,
data = bp, inflate = "bottom")
summary(m_ziop)
inflation_test(m_ziop) # Vuong (raw / AIC / BIC) vs the plain ordered probit
compare_models(op = m_op, ziop = m_ziop)
## regime-specific quantities
head(predict(m_ziop, type = "inflated")) # P(structurally peaceful)
predict(m_ziop, newdata = data.frame(loggdppc = 8, parliament = 1, disaster = 0,
major_oil = 0, major_primary = 0), se.fit = TRUE)
first_difference(m_ziop, "loggdppc", from = 7, to = 9) # total effect
first_difference(m_ziop, "loggdppc", from = 7, to = 9, stage = "inflation")
plot(ame(m_ziop, vars = c("loggdppc", "disaster"))) # average marginal effects
## correlated errors (ZiOPC)
m_ziopc <- iop(violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
loggdppc + parliament + disaster + major_oil + major_primary,
data = bp, inflate = "bottom", correlated = TRUE)
lr_test(m_ziop, m_ziopc) # rho = 0?
## top-inflated: escape-flexibility provisions in trade agreements
data(pta)
m_tiop <- iop(flexibility ~ depth * democracy + gdp + gdppc + trade + gattwto + members +
democratization | gdp + gdppc + democracy + democratization,
data = pta, inflate = "top")
## partial proportional odds for the plain ordered logit
m_pl <- ologit(flexibility ~ depth + democracy + gdp, data = pta, parallel = "auto")
m_pl$autofit$relaxed
## panels: cluster-robust SEs, random intercepts, the Mundlak device
m_cl <- iop(violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster,
data = bp, inflate = "bottom", cluster = "country")
m_re <- oprobit(violence ~ loggdppc + parliament + disaster, data = bp, re = "country")
head(ranef(m_re))
md <- mundlak(violence ~ loggdppc + disaster, data = bp, unit = "country")
oprobit(md$formula, md$data, cluster = "country")
## tables
texreg::screenreg(list(m_op, m_ziop))
modelsummary::modelsummary(list(OP = m_op, ZiOP = m_ziop))Outcome stage (ordered probit/logit with cutpoints tau):
P(y <= j | s = 1) = F(tau_j - x'beta). Inflation stage:
P(s = 1) = F(z'gamma) is the probability of the
ordered (non-inflated) regime, so positive inflation
coefficients raise the probability of the ordered regime; units in the
inflated regime (s = 0) are observed in the inflated
category k with certainty. Hence
P(y = j) = F(z'gamma) * pi_j(x) + 1{j = k} * [1 - F(z'gamma)].
With correlated = TRUE (probit only) the two latent
errors are bivariate normal with correlation rho, and the
regime-1 cell probabilities are rectangle probabilities of the bivariate
normal. See vignette("model").
The package’s tests check the bivariate-normal routine against
pbivnorm and mvtnorm, every analytic gradient
against numerical differentiation,
oprobit()/ologit() against
MASS::polr, ordinal::clm, and
VGAM::vglm, ordinal::clmm for the
random-intercept models, parameter recovery on simulated inflated data,
and – on the bundled bp data – Table 1 of Bagozzi, Hill,
Moore and Mukherjee (2015): OP, ZiOP, ZiOPC, and ZiOPC2 log-likelihoods
and the error correlations. Beyond the test suite, iop()
reproduces Table 1 and the marginal-effects figures of Bagozzi and
Mukherjee (2012), the TiOP/TiOPC log-likelihoods of Bagozzi, Joo and
Mukherjee (2024), and the Python idcempy implementation
where it converges; the original likelihood code of each paper evaluated
at iop’s estimates returns iop’s
log-likelihood to machine precision. See
data-raw/oracles/README.md in the source repository, and
inst/mc/ for the Monte Carlo behind the panel
recommendations.
citation("iop") gives the package citation and the
articles that introduced the models.
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.
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.