GLMM Inference for KK Designs with Ordinal Response
Source:R/inference_ordinal_KK_combined.R
InferenceOrdinalKKGLMM.RdFits a cumulative-logit random-intercept mixed model (proportional odds) for
ordinal responses under a KK matching-on-the-fly design:
\(\mathrm{logit}(P(Y_i \le k \mid w_i, x_i, b_{g(i)})) = \alpha_k - (\beta_T
w_i + x_i^\top \gamma) - b_{g(i)}\), for cutpoints \(\alpha_1 < \cdots <
\alpha_{K-1}\), treatment indicator \(w_i\), covariates \(x_i\), and a
matched-pair random intercept \(b_g \sim N(0, \sigma_b^2)\) that is
integrated out of the marginal likelihood (either by adaptive Gauss-Hermite
quadrature when use_rcpp = TRUE, the default; see
fast_ordinal_glmm_cpp for the quadrature order and optimizer
details, or by glmmTMB's Laplace approximation when use_rcpp =
FALSE). \(g(i)\) is subject \(i\)'s matched-pair group id; reservoir
(unmatched) subjects each get their own singleton group, contributing no
within-group correlation but still entering the joint likelihood. The
treatment coefficient \(\beta_T\) is a log-odds-ratio: \(\exp(\beta_T)\)
is the (conditional-on-\(b_g\)) odds ratio of being at or above any given
response category. likelihood_tier = "full": likelihood-ratio, score,
Wald, and gradient tests are available when use_rcpp = TRUE and the
fit converges; use_rcpp = FALSE disables likelihood-test support
(private$supports_likelihood_tests() returns FALSE) because
glmmTMB's Laplace-approximate likelihood is not wired into this
package's score/gradient/LR machinery. Validity requires the random-intercept
structure to correctly capture the design's matching dependence, proportional
odds (the treatment/covariate effect is constant across cutpoints), and
correct specification of the fixed-effects formula.
This differs from the GEE sibling
InferenceOrdinalKKGEE (documented
above) in estimand and inference basis: the GLMM's \(\beta_T\) is a
subject-specific (conditional) log-odds-ratio with model-likelihood-based
inference, while the GEE's is a population-averaged (marginal) log-odds-ratio
with sandwich-based inference; the two need not numerically agree even on the
same data, and the correct choice depends on whether a
subject-specific/conditional or population-averaged/marginal treatment
effect is of interest.
References
Hedeker, D., and Gibbons, R. D. (1994). "A Random-Effects Ordinal Regression Model for Multilevel Analysis." Biometrics, 50(4), 933-944, doi:10.2307/2533433 , for the random-effects cumulative-logit model; 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 the adaptive Gauss-Hermite quadrature approximation used to integrate out the random intercept.
See also
Comparable Python API: statsmodels MixedLM (continuous analog; no ordinal-GLMM in statsmodels). See also: Ordinal regression and Mixed model (Wikipedia).
Super class
Inference -> InferenceOrdinalKKGLMM
Methods
+ inherited public methods from Inference
Inference$capabilities()Inference$compute_exact_confidence_interval()Inference$compute_exact_two_sided_pval_for_treatment_effect()Inference$duplicate()Inference$get_analysis_data()Inference$get_covariates()Inference$get_design_object()Inference$get_model_formula()Inference$get_nonestimable_reason()Inference$get_nonestimable_stage()Inference$get_optimization_alg()Inference$get_response()Inference$get_response_type()Inference$get_treatment()Inference$is_nonestimable()Inference$set_optimization_alg()Inference$set_seed()Inference$supports()
InferenceOrdinalKKGLMM$new()
Initialize inference for the cumulative-logit random-intercept
mixed model \(\mathrm{logit}(P(Y_i \le k)) = \alpha_k - (\beta_T W_i +
X_i^\top \gamma) - b_{g(i)}\), \(b_g \sim N(0, \sigma_b^2)\), where
\(g(i)\) is subject \(i\)'s matched-pair group id (reservoir subjects
get singleton groups). Does not fit the model; the fit is deferred to
the first call to compute_estimate() or a method that requires it.
Usage
InferenceOrdinalKKGLMM$new(
des_obj,
model_formula = NULL,
use_rcpp = TRUE,
verbose = FALSE,
smart_cold_start_default = NULL
)Arguments
des_objA completed
Designobject with an ordinal response.model_formulaOptional formula for covariate adjustment. If
NULL(default), the formula from the design object is used and its pre-computed design matrix is reused. If a formula is provided, a new design matrix is constructed from the design's imputed covariates.use_rcppLogical. If
TRUE(default), use internal Rcpp.verboseWhether to print progress messages.
smart_cold_start_defaultWhether to use smart cold start values.
InferenceOrdinalKKGLMM$compute_estimate()
Fits the cumulative-logit random-intercept mixed model by
(adaptive-Gauss-Hermite- or Laplace-)approximate maximum likelihood and
returns \(\hat\beta_T\), the estimated treatment log-odds-ratio,
conditional on the matched-pair random intercept. Caches the fitted
model object, full parameter vector, and (when estimate_only =
FALSE) the standard error and degrees of freedom for reuse by
compute_asymp_confidence_interval(),
compute_asymp_two_sided_pval(), and likelihood-test methods; a
fit that fails the kernel's projected-gradient convergence check,
produces non-finite parameters, reaches the upper random-effect variance
boundary, exceeds private$max_abs_reasonable_coef, or lacks a
finite positive treatment-coefficient variance is cached as nonestimable
rather than returned. A valid near-zero random-effect variance boundary
is accepted using conditional fixed-effect information. The native
optimizer retains multistart L-BFGS for basin selection and, only when
its finite selected point fails the projected-score tolerance, applies
a damped-Newton polish using the numerical Hessian. The polished point
is retained only if it remains finite and does not increase the
negative log-likelihood; at a valid lower variance boundary, the
KKT-satisfied variance coordinate is excluded from that Newton system.
InferenceOrdinalKKGLMM$compute_asymp_confidence_interval()
Wald confidence interval for \(\beta_T\) using the fitted
model's standard error and degrees of freedom; see
InferenceAsymp for the shared
\(\hat\beta_T \pm t_{\alpha/2, df} \cdot \widehat{se}(\hat\beta_T)\)
(or z-based when df = Inf) contract. Fits the model first if not
already cached.
InferenceOrdinalKKGLMM$compute_asymp_two_sided_pval()
Two-sided Wald test of \(H_0: \beta_T = \code{delta}\)
against \(H_1: \beta_T \ne \code{delta}\), using the fitted model's
standard error and degrees of freedom; see
InferenceAsymp for the shared
\(t\)/\(z\) test contract. Fits the model first if not already
cached.
InferenceOrdinalKKGLMM$compute_estimate_with_bootstrap_weights()
Refits the mixed model with subject/block-level weights
applied to each row's contribution to the marginal likelihood
(Bayesian-bootstrap or nonparametric-bootstrap draw weights, expanded
from subject/block level to individual rows via
private$expand_subject_or_block_weights_to_row_weights()) and
returns the reweighted estimate \(\hat\beta_T^{(w)}\). Uses
fast_ordinal_regression_weighted_cpp — an ordinary
(non-mixed-effects) weighted cumulative-logit fit, not a reweighted
GLMM refit — as a fast approximation to the weighted marginal
likelihood; this trades exact random-effects refitting for speed across
many bootstrap replicates. When weights are effectively constant, this
collapses to the unweighted compute_estimate() call (returns
df = Inf to signal a degenerate/skipped bootstrap replicate
rather than refitting). Rows with non-finite or non-positive weight, or
non-finite response, are dropped from the weighted fit; if no rows
remain, the estimate is NA.
Examples
# \donttest{
seq_des = DesignSeqOneByOneKK14$new(n = 10, response_type = 'ordinal')
for (i in 1:10) {
seq_des$add_one_subject_to_experiment_and_assign(data.frame(x1 = rnorm(1), x2 = rnorm(1)))
}
seq_des$add_all_subject_responses(sample(1:4, 10, replace = TRUE))
inf = InferenceOrdinalKKGLMM$new(seq_des)
inf$compute_estimate()
#> [1] -0.3583555
# }