Fast Ordinal Cumulative-Logit Random-Intercept GLMM via Gauss-Hermite Quadrature (C++)

Fits a cumulative-logit ordinal mixed model with a single Gaussian random intercept per group (e.g. a matched pair or singleton from a KK-style matched design): $$\mathrm{logit}\,\Pr(Y_{ij} \le k \mid u_i) = \alpha_k - x_{ij}^\top\beta - u_i, \qquad u_i \sim N(0, \sigma^2),$$ for group \(i\), member \(j\), ordinal outcome \(y_{ij} \in \{1, \ldots, K\}\), and increasing cutpoints \(\alpha_1 < \cdots < \alpha_{K-1}\) (X carries no separate intercept column — the cutpoints serve that role). The marginal likelihood for group \(i\) integrates the random intercept out, $$L_i(\theta) = \int \prod_{j \in i} \Pr(Y_{ij} = y_{ij} \mid u_i) \, \phi(u_i / \sigma) \, du_i,$$ approximated by n_gh-point Gauss-Hermite quadrature (substituting \(u = \sqrt{2}\,\sigma\, z\) for quadrature node \(z\)), and optimized on the log scale by directly maximizing \(\sum_i \log L_i(\theta)\) (via optimize_fixed_likelihood, default optimization_alg = "lbfgs") over the reparameterized vector \([\alpha_1, \log(\alpha_2-\alpha_1), \ldots, \log(\alpha_{K-1}-\alpha_{K-2}), \beta, \log\sigma]\) — cutpoints are recovered as successive partial sums of \(\alpha_1\) and the exponentiated log-differences, which enforces \(\alpha_1 < \cdots < \alpha_{K-1}\) by construction rather than as a fitting constraint. Rows are stably sorted by group_id inside the kernel, so matched-group membership is invariant to input row order. Optimization uses supplied/cold, moderate-variance, and near-zero-variance starts, retains the smallest finite negative log-likelihood, and polishes that solution before applying a projected-gradient convergence check. If finite multistart L-BFGS stops on function decrease while its projected score remains above max(1e-5, eps_g), the kernel performs a local damped-Newton polish using its numerical Hessian. A Newton trial is retained only when its parameters, objective, and gradient are finite and its objective does not exceed the L-BFGS objective. At a valid lower log_sigma boundary, the KKT-satisfied variance coordinate is excluded from the Newton system, so fixed-effect convergence can be established without rejecting a near-zero random-effect variance. log_sigma is evaluated within \([-\code{max\_abs\_log\_sigma}, \code{max\_abs\_log\_sigma}]\) with a quadratic penalty on excursions beyond that interval whose analytic gradient matches the bounded objective; variance_boundary_hit in the return value flags whether the fitted log_sigma landed at that boundary (a sign the random-intercept variance is being driven to (near-)zero or is unbounded, and that ssq_b_T/fisher_information should be treated with caution). The Hessian used for inference is a numerical (central finite-difference, step \(10^{-4}\), symmetrized) second derivative of the analytic gradient, not a closed-form expression. At the valid near-zero variance boundary, treatment variance is computed conditional on that boundary by excluding the nonregular variance-parameter row and column.

Usage

fast_ordinal_glmm_cpp(
  X,
  y,
  group_id,
  K,
  j_T,
  smart_cold_start = TRUE,
  estimate_only = FALSE,
  n_gh = 20L,
  max_abs_log_sigma = 8,
  maxit = 300L,
  eps_g = 1e-06,
  warm_start_params = NULL,
  warm_start_beta = NULL,
  optimization_alg = "lbfgs",
  fixed_idx = NULL,
  fixed_values = NULL,
  warm_start_fisher_info = NULL
)

Arguments

X

A numeric matrix of predictors, one row per observation (member-level, not group-level); no intercept column (see Details).

y

Integer vector of 1-indexed ordinal outcomes (\(1, \ldots, K\)), one per row of X.

group_id

Integer vector of group (e.g. matched-pair) identifiers, one per row of X; assumed contiguous per group after internal sorting-by-value into blocks.

K

The number of ordinal levels.

j_T

0-based column index of X whose coefficient's variance (ssq_b_T) should be computed — typically the treatment indicator.

smart_cold_start

Logical. If TRUE and no warm start is given, initialize cutpoints evenly at 0 (all log-differences zero) and \(\beta\) via a naive OLS fit of y (treated as numeric) on X; log_sigma starts at \(-3\). If FALSE, \(\beta\) instead starts at zero.

estimate_only

If TRUE, skip the (converged-fit-only) covariance calculation for ssq_b_T — point estimates and the Hessian are still returned regardless.

n_gh

Number of Gauss-Hermite quadrature nodes used to integrate out the random intercept.

max_abs_log_sigma

Symmetric clamp bound for log_sigma during optimization (default 8).

maxit

Maximum number of optimizer iterations.

eps_g

Gradient-norm convergence tolerance.

warm_start_params

Optional starting values for the full reparameterized vector \([\alpha_1, \log\text{-diffs}, \beta, \log\sigma]\); if its length doesn't match, falls back to zero cutpoints/\(\beta\) and log_sigma = -3. Takes precedence over warm_start_beta and smart_cold_start.

warm_start_beta

Optional starting values either for the full parameter vector (same length as warm_start_params above) or for \(\beta\) alone (length p, with cutpoints zeroed and log_sigma = -3); ignored if warm_start_params is supplied.

optimization_alg

Optimization algorithm (default "lbfgs").

fixed_idx

Optional 1-indexed positions (into the reparameterized parameter vector) to hold fixed.

fixed_values

Optional values, parallel to fixed_idx, of the fixed parameters.

warm_start_fisher_info

Optional initial curvature matrix for the first optimizer iteration.

Value

A list with components b (\(\hat\beta\)), alpha (the \(K-1\) cutpoints, recovered from the reparameterization), params (the full fitted reparameterized vector), log_sigma, ssq_b_T (variance of b[j_T], NA unless estimate_only = FALSE and the fit converged and the resulting information matrix inverts successfully), converged, neg_loglik, fisher_information (the numerical Hessian, always returned), score (the log-likelihood score at the returned parameters), gradient_norm, newton_polish_attempted, newton_polish_accepted, and newton_polish_iterations (diagnostics for the conditional damped-Newton fallback), and variance_boundary_hit (TRUE/FALSE, or NA if the optimizer itself threw an exception, in which case converged = FALSE and all other quantities besides b/alpha/log_sigma are NA).

References

Pinheiro, J. C., and Bates, D. M. (1995). "Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Effects Model." Journal of Computational and Graphical Statistics, 4(1), 12-35, doi:10.1080/10618600.1995.10474663 , for Gauss-Hermite quadrature as an approximation to the random-effect marginal likelihood integral used throughout this package's GLMM backends (fast_poisson_glmm_cpp, fast_logistic_glmm_cpp, fast_weibull_frailty_cpp, and this function).