Ordered probit and ordered logit regression

Fits the standard ordered probit (oprobit()) or ordered logit (ologit()) model by maximum likelihood, with optional partial proportional-odds (non-parallel, "generalized ordered") effects, a unit random intercept, survey weights, offsets, and analytic, robust, or cluster-robust standard errors. These are the matched baselines for the inflated models iop() and iol(): the same engine, the same cutpoint conventions, and the same methods, so comparisons across the four are like for like.

Usage

oprobit(
  formula,
  data,
  parallel = TRUE,
  re = NULL,
  nAGQ = 15,
  fe = NULL,
  fe_correction = c("none", "jackknife"),
  time = NULL,
  weights = NULL,
  offset = NULL,
  se = c("analytic", "robust", "cluster", "bootstrap", "none"),
  cluster = NULL,
  nboot = 200,
  cores = 1,
  start = NULL,
  maxit = 1000,
  reltol = 1e-10
)

ologit(
  formula,
  data,
  parallel = TRUE,
  re = NULL,
  nAGQ = 15,
  fe = NULL,
  fe_correction = c("none", "jackknife"),
  time = NULL,
  weights = NULL,
  offset = 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.
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.
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_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.
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("oprobit", "iord") or c("ologit", "iord"): a list with coefficients (outcome coefficients, then the cutpoints named "a|b", then sigma_u for a random-intercept fit), vcov, loglik, fitted.values (the n x J matrix of category probabilities), cutpoints, converged, boundary, and the design pieces used by the methods. See summary.iord(), predict.iord(), first_difference(), compare_models().

Details

The latent-variable model is $y_i^* = x_i'\beta + \epsilon_i$ with $y_i = j$ when $\tau_{j-1} < y_i^* \le \tau_j$ ($\tau_0 = -\infty$, $\tau_J = \infty$), so $P(y_i \le j) = F(\tau_j - x_i'\beta)$, with $F$ the standard normal (probit) or logistic (logit) distribution function. The outcome equation has no intercept (it is absorbed by the cutpoints), matching MASS::polr().

With parallel = FALSE, or a one-sided formula naming the terms held parallel (all others relaxed; ~ . - x2 relaxes only x2), the relaxed terms get 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). Such fits can imply negative cell probabilities when the category-specific curves cross; the fit warns when that happens at the optimum.

With re = "unit" a random intercept $u_g \sim N(0, \sigma_u^2)$ for the units of that column enters the latent equation, $y_{ig}^* = x_{ig}'\beta + u_g + \epsilon_{ig}$, and is integrated out by adaptive Gauss–Hermite quadrature with nAGQ nodes per unit (the mode and curvature of each unit's integrand re-center and re-scale the nodes), as in ordinal::clmm(). Reported probabilities and first differences are then population-averaged (marginal over the random intercept); conditional probabilities at $u = 0$ and the empirical-Bayes unit effects are available through predict.iord() and ranef().

For panels, mundlak() adds unit means of the covariates (correlated random effects), which pairs naturally with re = or with cluster =.

Missing values, refits, and new data

Rows with a missing value in any variable of either equation (or in weights, offsets, cluster, re, fe) are dropped before fitting – there is no na.action argument – and object$data, fitted(), residuals(), and predict() without newdata refer to the retained rows (so they are not padded to the original row positions as with na.exclude). The stored call makes update.iord() work as usual (update(fit, . ~ . + x3), update(fit, data = subset); for the two-part formulas update(fit, . ~ . + x3 | .) changes the outcome equation and update(fit, . ~ . | . + z2) the inflation equation). In predict.iord(), newdata must contain every covariate of both equations; extra columns are ignored, unseen factor levels are an error, and rows with missing covariates give NA predictions.

References

McKelvey, R.D. and Zavoina, W. (1975). A statistical model for the analysis of ordinal level dependent variables. Journal of Mathematical Sociology, 4, 103-120.

Peterson, B. and Harrell, F.E. (1990). Partial proportional odds models for ordinal response variables. Applied Statistics, 39, 205-217.

Williams, R. (2006). Generalized ordered logit/partial proportional odds models for ordinal dependent variables. Stata Journal, 6, 58-82.

Greene, W. (2004). The behaviour of the maximum likelihood estimator of limited dependent variable models in the presence of fixed effects. Econometrics Journal, 7, 98-119.

Dhaene, G. and Jochmans, K. (2015). Split-panel jackknife estimation of fixed-effect models. Review of Economic Studies, 82, 991-1030.

Dale, D. and Sirchenko, A. (2021). Estimation of nested and zero-inflated ordered probit models. Stata Journal, 21, 3-38.

Examples

set.seed(1)
n <- 500; x <- rnorm(n); z <- rbinom(n, 1, 0.5)
ystar <- 0.8 * x - 0.5 * z + rnorm(n)
y <- cut(ystar, c(-Inf, -0.7, 0.5, Inf), labels = c("low", "mid", "high"), ordered_result = TRUE)
d <- data.frame(y, x, z)
m1 <- oprobit(y ~ x + z, data = d)
summary(m1)
#> 
#> Ordered probit
#> Call:  oprobit(formula = y ~ x + z, data = d)
#> Response levels (in order): low < mid < high 
#> N = 500   inference: analytic 
#> 
#> Outcome equation (ordered probit):
#>   Estimate Std. Error z value Pr(>|z|)    
#> x  0.79031    0.06294  12.556  < 2e-16 ***
#> z -0.41551    0.10750  -3.865 0.000111 ***
#> 
#> Cutpoints:
#>          Estimate Std. Error z value Pr(>|z|)    
#> low|mid  -0.67575    0.08277  -8.164 3.25e-16 ***
#> mid|high  0.50641    0.08116   6.239 4.39e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> logLik = -453.26   AIC = 914.52   BIC = 931.37   df = 4
m2 <- ologit(y ~ x + z, data = d, se = "robust")
compare_models(probit = m1, logit = m2)
#>    model           type    logLik df      AIC      BIC   N
#> 1 probit Ordered probit -453.2582  4 914.5164 931.3748 500
#> 2  logit  Ordered logit -453.6283  4 915.2565 932.1150 500
head(predict(m1, type = "prob"))
#>             low       mid       high
#> [1,] 0.59283812 0.3289222 0.07823965
#> [2,] 0.34260009 0.4387566 0.21864332
#> [3,] 0.49387755 0.3844796 0.12164289
#> [4,] 0.06412858 0.3032328 0.63263867
#> [5,] 0.17459334 0.4225627 0.40284399
#> [6,] 0.48909882 0.3868216 0.12407962