---
title: "Getting started with iop"
author: "Benjamin E. Bagozzi"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Getting started with iop}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4,
                      message = FALSE)
```

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

```{r bp}
library(iop)
data(bp)
round(prop.table(table(bp$violence)), 3)
```

The 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):

```{r bp-fit}
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)
```

The 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:

```{r bp-test}
inflation_test(m_ziop)
compare_models(op = m_op, ziop = m_ziop)
```

### Quantities 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:

```{r bp-pred}
head(predict(m_ziop, type = "prob"), 3)
summary(predict(m_ziop, type = "inflated"))
summary(predict(m_ziop, type = "posterior")[bp$violence == "none"])
```

First 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:

```{r bp-fd}
first_difference(m_ziop, "loggdppc", from = 7, to = 9)
first_difference(m_ziop, "loggdppc", from = 7, to = 9, stage = "inflation")
```

`vignette("quantities")` covers predictions on covariate profiles, standard
errors for predicted probabilities, average marginal effects, plots, and
regression tables.

### Correlated errors

`correlated = TRUE` estimates the correlation between the two latent
equations' errors (the ZiOPC model). The correlated fit starts from the
uncorrelated solution and nests it, so a likelihood-ratio test of `rho = 0` is
standard:

```{r bp-corr}
m_ziopc <- iop(f_ziop, data = bp, inflate = "bottom", correlated = TRUE)
coef(m_ziopc)["rho"]
confint(m_ziopc, parm = "rho")
lr_test(m_ziop, m_ziopc)
```

## 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:

```{r bp-logit}
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)
```

## Top 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:

```{r pta}
data(pta)
table(pta$flexibility)
m_tiop <- iop(flexibility ~ depth * democracy + gdp + gdppc + trade + gattwto + members +
                democratization | gdp + gdppc + democracy + democratization,
              data = pta, inflate = "top")
summary(m_tiop)
```

Here 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:

```{r repression}
data(repression)
m_rep <- iop(repression ~ negxpol * oilrent | dom_media + civil_war,
             data = repression, inflate = "top")
inflation_test(m_rep)$vuong
```

## Middle 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:

```{r middle}
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)
m_miop <- iop(y ~ x1 + x2 | z1, data = d, inflate = "middle")
m_miop
```

The 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")`:

```{r middle-post}
summary(predict(m_miop, type = "posterior")[d$y == "neutral"])
```

`inflate` 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):

```{r affairs, eval = requireNamespace("AER", quietly = TRUE)}
data("Affairs", package = "AER")
m_aff <- iop(affairs ~ age + yearsmarried + religiousness + rating | gender + children,
             data = Affairs, inflate = "bottom")
m_aff
colMeans(predict(m_aff, type = "zeros"))       # zeros from the inflation process / the ordered stage
inflation_test(m_aff)
```

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:

```{r nmes, eval = requireNamespace("AER", quietly = TRUE)}
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)
ame(m_nmes, vars = c("chronic", "insurance"))
```

Everything 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 of `split = "category"` and
  `split_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.