---
title: "The model, identification, and estimation"
author: "Benjamin E. Bagozzi"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{The model, identification, and estimation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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.

## The ordered stage

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.

### Partial proportional odds

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

```{r parallel}
library(iop)
data(pta)
m_pl <- ologit(flexibility ~ depth + democracy + gdp + gdppc + trade, data = pta)
parallel_test(m_pl)
m_auto <- ologit(flexibility ~ depth + democracy + gdp + gdppc + trade, data = pta, parallel = "auto")
m_auto$autofit$relaxed
lr_test(m_pl, m_auto)
```

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

## The inflation stage

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")`).

### Correlated errors

With `correlated = TRUE` (probit only), \((u_i, \epsilon_i)\) are bivariate
standard normal with correlation \(\rho\), and the regime-1 cell
probabilities become rectangle probabilities of the bivariate normal,
\(P(s_i = 1, y_i \le j) = \Phi_2(z_i'\gamma, \tau_j - x_i'\beta; -\rho)\),
evaluated by a deterministic Gauss--Legendre algorithm (Drezner and Wesolowsky
1990; Genz 2004) so that the likelihood is smooth for the optimizer. \(\rho\)
is estimated on the \(\tanh^{-1}\) scale and reported on the natural scale;
`confint()` transforms the interval back so it respects \((-1, 1)\). The
correlated fit starts from the uncorrelated solution and nests it, so
`lr_test()` of \(\rho = 0\) is a standard one-degree-of-freedom test:

```{r tiopc}
f <- flexibility ~ depth * democracy + gdp + gdppc + trade + gattwto + members + democratization |
  gdp + gdppc + democracy + democratization
m_tiop  <- iop(f, data = pta, inflate = "top")
m_tiopc <- iop(f, data = pta, inflate = "top", correlated = TRUE)
c(rho = coef(m_tiopc)["rho"], confint(m_tiopc, parm = "rho"))
lr_test(m_tiop, m_tiopc)
```

There is no canonical bivariate logistic distribution, so `iol()` has no
correlated option.

### Category-specific split equations

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

```{r gziop, message = FALSE}
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")
compare_models(common = m_common, category = m_category)
split_test(m_common)
m_category
```

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"`.

## Identification

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:

```{r exclusion, message = TRUE}
data(bp)
m_same <- iop(violence ~ loggdppc + disaster | loggdppc + disaster, data = bp, inflate = "bottom")
```

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.

## Estimation

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.

### Multi-start and the ordered baseline

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:

```{r multistart, message = FALSE}
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)
```

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

### Boundary cases

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:

```{r boundary, warning = TRUE}
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
m0 <- iop(y ~ x1 + x2 | z1, data = d0, inflate = "bottom")
m0$boundary
c(inflated = m0$loglik, ordered = m0$loglik_uninflated)
```

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

## Which test for which comparison

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

```{r boot-lr, eval = FALSE}
inflation_test(m_ziop, boot = 199)
```

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, offsets, standard errors

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

## Names in the literature

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 |

## Validation

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.

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

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.
