The model, identification, and estimation

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:

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.

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:

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

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:

data(bp)
m_same <- iop(violence ~ loggdppc + disaster | loggdppc + disaster, data = bp, inflate = "bottom")
#> 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.

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:

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

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.

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:

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.

Benjamin E. Bagozzi