Inflated ordered probit and inflated ordered logit regression

Fits an ordered probit (iop()) or ordered logit (iol()) model in which one ordered category is inflated: it collects observations generated by the ordered process together with observations generated by a distinct split-population process that places them in that category regardless of the ordered mechanism. The zero-inflated ordered probit of Harris and Zhao (2007), the middle-inflated ordered probit of Bagozzi and Mukherjee (2012), and the top-inflated ordered probit of Bagozzi, Joo and Mukherjee (2024) are the cases inflate = "bottom", "middle", and "top"; any single category may be named. iop() optionally estimates the correlation between the two latent equations' errors (correlated = TRUE, the ZiOPC / MiOPC / TiOPC models). Both accept unit random intercepts (re), in the outcome equation and optionally also in the inflation equation.

Usage

iop(
  formula,
  data,
  inflate,
  correlated = FALSE,
  split = c("common", "category"),
  parallel = TRUE,
  re = NULL,
  re_inflation = FALSE,
  nAGQ = 15,
  fe = NULL,
  fe_inflation = FALSE,
  fe_correction = c("none", "jackknife"),
  time = NULL,
  weights = NULL,
  offset = NULL,
  offset_inflation = NULL,
  se = c("analytic", "robust", "cluster", "bootstrap", "none"),
  cluster = NULL,
  nboot = 200,
  cores = 1,
  start = NULL,
  maxit = 1000,
  reltol = 1e-10
)

iol(
  formula,
  data,
  inflate,
  split = c("common", "category"),
  parallel = TRUE,
  re = NULL,
  re_inflation = FALSE,
  nAGQ = 15,
  fe = NULL,
  fe_inflation = FALSE,
  fe_correction = c("none", "jackknife"),
  time = NULL,
  weights = NULL,
  offset = NULL,
  offset_inflation = NULL,
  se = c("analytic", "robust", "cluster", "bootstrap", "none"),
  cluster = NULL,
  nboot = 200,
  cores = 1,
  start = NULL,
  maxit = 1000,
  reltol = 1e-10
)

Arguments

formula: A model formula y ~ x1 + x2. The response may be an ordered factor, a factor (level order taken as the ordinal order), or a numeric/integer vector (sorted unique values define the order). The outcome equation has no intercept (the cutpoints absorb it), so an explicit 0 + / - 1 is ignored with a message and factors are coded as with an intercept; a set of columns that sums to a constant is rejected as rank deficient.
data: A data frame.
inflate: The inflated category: "bottom", "middle" (odd number of categories only), "top", or one of the response's categories (a level label or, for a numeric response, a value).
correlated: Logical (iop() only): estimate the correlation between the inflation- and outcome-equation errors.
split: "common" (default): one split equation shared by every category (the ZiOP / MiOP / TiOP); "category": a separate split equation for every non-inflated category (the generalised GZiOP / GMiOP of Brown, Harris and Spencer 2020; see Details).
parallel: TRUE (default; all effects parallel), FALSE (all effects category-specific), a one-sided formula naming the terms that are held parallel, e.g. parallel = ~ . - x2, or "auto": starting from the all-parallel fit, relax the term whose parallel restriction has the smallest likelihood-ratio p-value while that p-value is below 0.05 (a forward version of Stata's gologit2, autofit; the steps are stored in object$autofit). See also parallel_test().
re: Optional column name in data identifying units that receive a random intercept in the outcome equation.
re_inflation: Logical: with re, also give the inflation equation an (independent) unit random intercept.
nAGQ: Number of adaptive Gauss–Hermite quadrature nodes per random intercept (default 15).
fe: Optional column name in data identifying units that receive a fixed effect (a dummy per unit, reference level omitted) in the outcome equation. Units whose response is constant at an extreme category have no finite fixed effect and are dropped with a message; a covariate that does not vary within units is collinear with the dummies and is refused by name (keep it with re = or mundlak() instead). Maximum-likelihood unit dummies carry incidental-parameters bias when units have few observations (Greene 2004); the fit warns when the median unit has fewer than 10. In the package's Monte Carlo (system.file("mc", package = "iop"): ordered probit, 100 units, a covariate correlated with a N(0,1) unit effect, 100 replications) the bias of that covariate's coefficient is +26 / +11 / +5 / +2.5 percent at T = 4 / 8 / 16 / 32 with unit dummies, +21 / +11 / +6 / +3 percent with a random intercept, +12 percent throughout when pooled, and within 0.3 percent at every T with the mundlak() device – the recommended route for short panels; reserve fe = for long ones. Computationally the dummies enter the parameter vector one per unit, and the exact-Hessian Newton polish and the covariance cost grow roughly with the square of the number of parameters (about 0.5 s at 20 units, 3 s at 100, 18 s at 300 for a plain ordered probit on 4,000 rows; the multistart of an inflated model multiplies this); for panels with many hundreds of units prefer re = or mundlak().
fe_inflation: Logical: with fe, also give the inflation equation unit dummies; units never or always in the inflated category are then dropped. A unit whose inflated-category observations are all absorbed by the ordered stage still has an inflation dummy at the boundary (regime probability 1); the fit reports that dummy's standard error as NA ($se_na) and flags the conditioning ($ill_conditioned); like the outcome-equation dummies, the inflation dummies are exempt from the quasi-separation signatures.
fe_correction: With fe, "jackknife" applies the split-panel jackknife of Dhaene and Jochmans (2015): the model is refit on the two half-panels of every unit (first and second half of the observations ordered by time, or by row order if time is missing) and the common parameters are bias-corrected as 2 * full - (half1 + half2) / 2, removing the leading incidental-parameters bias; the uncorrected estimates are kept in object$coefficients_uncorrected, the half-panel estimates and log-likelihoods in object$jackknife, and the full-sample covariance is reported: Dhaene and Jochmans (2015, Section 2) show that the split-panel jackknife removes the leading bias without changing the first-order asymptotic variance, so the full-sample covariance is the asymptotically valid one for the corrected estimator (its finite-sample variance is somewhat larger, which se = "bootstrap" cannot assess under fe; treat the intervals as approximate). Units with a single observation do not enter the half-panels. The correction is only as good as the half-panel fits: when a covariate is identified mainly by a within-unit trend, or the half-panels are very short, the half-panel estimates are noisy and the corrected estimate inherits that noise – inspect object$jackknife$half_coefficients before reporting it.
time: Optional column name ordering the observations within units (used by fe_correction = "jackknife").
weights: Optional weights: a column name in data or a numeric vector aligned to its rows. They enter 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 not frequency weights: 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.
offset: Optional offset on the latent scale of the outcome equation: a column name in data or a numeric vector.
offset_inflation: Optional offset on the latent scale of the inflation equation: a column name in data or a numeric vector.
se: Standard errors: "analytic" (inverse observed information, default), "robust" (heteroskedasticity-consistent sandwich over observations, or over units for random-intercept fits), "cluster" (cluster-robust; needs cluster), "bootstrap" (nonparametric bootstrap over the estimation rows – or over the clusters when cluster is given, or over the units of a random-intercept fit – with nboot refits from the full-sample estimate; the replicate estimates are kept in object$boot and confint.iord() can then report percentile intervals), or "none". The bootstrap is not available with fe. Replicates whose split equation (or rho) runs to a boundary are counted in object$boot$n_boundary and trigger a warning: their estimates inflate the bootstrap standard errors of the affected block, for which the percentile intervals are the more robust summary. For "bootstrap", a seed set before the call (set.seed()) makes the resamples reproducible; runs with cores = 1 and cores > 1 use different random streams.
cluster: Optional cluster identifier (a column name in data or a vector aligned to its rows); supplying it selects se = "cluster" unless se is given explicitly (with se = "bootstrap" it defines the resampling blocks). With a random intercept, clusters must nest the units.
nboot: Number of bootstrap refits for se = "bootstrap" (default 200). Each refit costs about one fit of the model (roughly a third of a second per refit for the bundled bp ZiOP), so budget a minute or two at the default on data of that size, and use cores.
cores: Number of parallel workers for the bootstrap refits (a PSOCK cluster via the parallel package; default 1).
start: Optional starting values on the internal parameter scale (see object$theta); rarely needed.
maxit, reltol: Controls for the BFGS stage of the optimizer, passed to stats::optim(); the exact-Hessian Newton polish that follows runs up to 25 further iterations regardless of maxit (so maxit bounds the quasi-Newton phase, not the total). Results are insensitive to maxit above a few hundred; maxit = 0 is reported as non-converged.

Value

An object of class c("iop", "iord") or c("iol", "iord"), as for oprobit(), with additional components: inflate (the inflated category label), k (its 0-based index), regime (the fitted $P(s_i = 1)$), rho (when correlated = TRUE), sigma_u / sigma_v (random-intercept SDs), loglik_uninflated (the plain ordered baseline), boundary, and coefficients named infl_<term> for the inflation equation. The coefficient vector is ordered outcome terms, cutpoints ("a|b"), inflation terms, rho, sigma_u, sigma_v.

Details

Two latent equations are estimated jointly. The inflation (split) equation $s_i^* = z_i'\gamma + u_i$ assigns unit $i$ to the ordered regime ($s_i = 1$, 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. The outcome equation is the ordered model of oprobit()/ologit() with cell probabilities $\pi_{ij}$. Hence $$P(y_i = j) = F(z_i'\gamma)\,\pi_{ij} + 1\{j = k\}\,[1 - F(z_i'\gamma)].$$ Note the sign convention, which follows 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)$.

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 the likelihood is smooth for the optimizer. $\rho$ is estimated on the $\tanh^{-1}$ scale and reported on the natural scale; its confidence interval is transformed from the former.

The inflation equation is specified after a | in the formula, y ~ x1 + x2 | z1 + z2, and includes an intercept. If the | part is omitted the outcome covariates are reused, which identifies the model by functional form alone; an exclusion restriction (a covariate in one equation only) is advisable and a message says so.

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 – so every inflated fit is 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, and the best optimum is kept. Starts are staged: each gets a short quasi-Newton run, and the two best are run to convergence and polished with exact-Hessian Newton steps. object$start_logliks records the log-likelihood reached from every start (short-run values for the non-finalists); a spread of several units across starts is the signature of multimodality, and start = lets you add 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. Fits whose inflation equation degenerates – every regime probability at 1 (the fit has collapsed to the plain ordered model) or a quasi-separated split with a few units at 0 and the rest at 1 (the coefficients are not finite) – whose rho sits at $\pm 1$, or whose random-intercept SD is zero are flagged (object$boundary) with a warning that names the case; such data contain no identifiable inflation (or correlation) process, or an over-rich split equation.

Random intercepts (re) enter the outcome equation; with re_inflation = TRUE an independent random intercept also enters the inflation equation, integrated by two-dimensional adaptive Gauss–Hermite quadrature (nAGQ^2 nodes per unit).

Diagnostics and fixed thresholds

Every fit carries programmatic diagnostics that summary() and print() also report. object$boundary (with object$boundary_messages) is TRUE when (i) a split equation is degenerate – every fitted regime probability within 1e-6 of 0 or 1 (collapse to the plain ordered model, or a perfect split); (ii) a split equation is quasi-separated: a standardized split slope (coefficient times the covariate's SD) above 10 or a split intercept above 40 in absolute value, i.e. the split acts as a step function and its maximum-likelihood estimate may not be finite (legitimately sharp splits in the package's applications stay below 6 and 21), or a split standard error above 50 on the standardized scale or below 1e-6 times its coefficient; the same signatures (standardized coefficient above 10, or standardized standard error above 50) flag a quasi-separated outcome equation, unit dummies excepted; (iii) rho beyond 0.985 in absolute value; (iv) a random-intercept SD below 1e-3; (v) the ordered stage supplies less than 1e-6 of the inflated category's fitted probability mass, or less than 1e-2 of it while the adjacent cutpoint's standard error exceeds 50 – the inflation process absorbs every observation in that category (a hurdle-type split), so the adjacent cutpoint is at $-\infty$ (bottom) or $+\infty$ (top), or the two cutpoints bracketing a middle category coincide; the finite value the optimizer stops at is an artifact and that cutpoint's standard error is reported as NA. Degenerate split- or outcome-equation standard errors are left visible (they are the evidence of the problem and the warning names them); only a cutpoint at $\pm\infty$ is reported as NA, because the finite value the optimizer stops at carries no information. The flags also travel with glance() and the texreg GOF block. object$converged is FALSE when the Newton decrement at the optimum exceeds 1e-6 and the gradient (per observation) exceeds 1e-6 or the Hessian is indefinite; object$max_grad is the largest gradient component at the optimum in the optimizer's standardized parameterization (the units the convergence rule uses; the raw gradient scales with the covariates' units). converged is a statement about the optimizer – the returned point is a stationary point of the likelihood – while the boundary flags are statements about that point; a quasi-separated fit can be converged in this sense, so check both. object$flat_hessian marks a converged optimum with a near-flat direction (relative eigenvalue of the final scaled Hessian below 1e-10, weak identification; every ill-conditioned fit is flat, a flat fit whose information matrix still inverts has computable but large standard errors); object$ill_conditioned marks an information matrix with reciprocal condition number below 1e-12 (near-collinear covariates or a near-flat split), computed in the optimizer's scaled parameterization – every covariate column standardized – so that the units of the covariates do not affect it (the Hessian itself is differenced in that parameterization, so the standard errors, like the condition number, are invariant to a change of units); object$se_na lists parameters whose delta-method variance is not positive (reported as NA, never as 0). These thresholds are fixed by design (they are not tuning parameters) and were chosen from the package's applications and stress tests so that every legitimate fit in them passes and every constructed boundary case is caught; see vignette("model").

Category-specific split equations. With split = "category" every non-inflated category $j$ gets its own split equation $s_{ij}^* = z_i'\gamma_j + u_{ij}$ (and, with correlated = TRUE, its own $\rho_j$): a unit whose ordered outcome would be $j$ is "tempered" into the inflated category with probability $1 - F(z_i'\gamma_j)$, so $$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}.$$ This is the generalised zero-/middle-inflated ordered probit (GZiOP, GMiOP, and their correlated versions) of Brown, Harris and Spencer (2020), which nests the common-split model when all $\gamma_j$ (and $\rho_j$) are equal; split_test() tests that restriction by a score (LM) test from the common-split fit and by the likelihood-ratio test. The fit starts from its common-split counterpart (kept in object$loglik_common) and checks that it is not below it. Identification now rests on $J - 1$ split equations, so exclusion restrictions matter more; coefficients are named infl_<term>:<category> and rho:<category>, and the regime quantities (predict(type = "regime"), $regime) become one column per non-inflated category. Not available with random intercepts or fe_inflation.

References

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.

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. Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing, 14, 251-260.

Examples

set.seed(2)
d <- riop(n = 800, beta = c(0.8, -0.5), gamma = c(0.4, 1), tau = c(-0.6, 0.7),
          inflate = "bottom")
m <- iop(y ~ x1 + x2 | z1, data = d, inflate = "bottom")
summary(m)
#> 
#> Inflated ordered probit (inflated category: 0)
#> Call:  iop(formula = y ~ x1 + x2 | z1, data = d, inflate = "bottom")
#> Response levels (in order): 0 < 1 < 2 
#> N = 800   inference: analytic 
#> 
#> Outcome equation (ordered probit):
#>    Estimate Std. Error z value Pr(>|z|)    
#> x1  0.74425    0.06884  10.812  < 2e-16 ***
#> x2 -0.45010    0.06767  -6.651 2.91e-11 ***
#> 
#> Cutpoints:
#>     Estimate Std. Error z value Pr(>|z|)    
#> 0|1  -0.5946     0.1161  -5.121 3.03e-07 ***
#> 1|2   0.6518     0.0709   9.192  < 2e-16 ***
#> 
#> Inflation equation (P(ordered regime); inflated category "0", observed share 0.562, mean fitted P(ordered regime) 0.641):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   0.4614     0.1234   3.739 0.000185 ***
#> z1            0.9602     0.1158   8.295  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> logLik = -635.41   AIC = 1282.82   BIC = 1310.93   df = 6
head(predict(m, type = "prob"))
#>              0         1          2
#> [1,] 0.5425901 0.3397206 0.11768939
#> [2,] 0.4417293 0.1544114 0.40385927
#> [3,] 0.3108475 0.2182638 0.47088874
#> [4,] 0.5808083 0.3418365 0.07735523
#> [5,] 0.6905407 0.2034067 0.10605257
#> [6,] 0.3786220 0.3226660 0.29871199
head(predict(m, type = "inflated"))         # P(inflated regime)
#> [1] 0.21572706 0.41873409 0.27302880 0.04765278 0.56322269 0.26280245
m0 <- oprobit(y ~ x1 + x2, data = d)
vuong(m, m0)                                # inflated vs plain ordered probit
#> Vuong test:  m1 = m [ Inflated ordered probit (inflated category: 0) ]
#>              m2 = m0 [ Ordered probit ]
#> positive statistics favor m1
#> 
#>  correction statistic p_m1_better p_m2_better p_two_sided
#>         raw     7.492           0           1           0
#>         AIC     7.334           0           1           0
#>         BIC     6.965           0           1           0
# \donttest{
mc <- iop(y ~ x1 + x2 | z1, data = d, inflate = "bottom", correlated = TRUE)
summary(mc)
#> 
#> Inflated ordered probit (inflated category: 0), correlated errors
#> Call:  iop(formula = y ~ x1 + x2 | z1, data = d, inflate = "bottom", 
#>     correlated = TRUE)
#> Response levels (in order): 0 < 1 < 2 
#> N = 800   inference: analytic 
#> 
#> Outcome equation (ordered probit):
#>    Estimate Std. Error z value Pr(>|z|)    
#> x1  0.73483    0.07025  10.461  < 2e-16 ***
#> x2 -0.43971    0.06808  -6.458 1.06e-10 ***
#> 
#> Cutpoints:
#>     Estimate Std. Error z value Pr(>|z|)    
#> 0|1  -0.6527     0.1355  -4.818 1.45e-06 ***
#> 1|2   0.5599     0.1220   4.591 4.41e-06 ***
#> 
#> Inflation equation (P(ordered regime); inflated category "0", observed share 0.562, mean fitted P(ordered regime) 0.650):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   0.4883     0.1302   3.751 0.000176 ***
#> z1            0.9197     0.1226   7.503 6.22e-14 ***
#> 
#> Error correlation:
#>     Estimate Std. Error z value Pr(>|z|)
#> rho  -0.2359     0.2307  -1.022    0.307
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> logLik = -634.91   AIC = 1283.81   BIC = 1316.60   df = 7
# }