The model, identification, and estimation

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:

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

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.

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:

f <- flexibility ~ depth * democracy + gdp + gdppc + trade + gattwto + members + democratization |
  gdp + gdppc + democracy + democratization
m_tiop  <- iop(f, data = pta, inflate = "top")
#> 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.
m_tiopc <- iop(f, data = pta, inflate = "top", correlated = TRUE)
#> 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.
c(rho = coef(m_tiopc)["rho"], confint(m_tiopc, parm = "rho"))
#>     rho.rho                         
#> -0.01486501 -0.92973338  0.92558846
lr_test(m_tiop, m_tiopc)
#> Likelihood-ratio test
#>   LR = 0.0003 on 1 df, p = 0.9861

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

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

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:

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

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:

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

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.