Fits ordinal cumulative probability models for continuous or ordinal
response variables, efficiently allowing for a large number of
intercepts by capitalizing on the information matrix being sparse.
Five different distribution functions are implemented, with the
default being the logistic (yielding the proportional odds model).
Penalized estimation and weights are also implemented, as in
lrm.fit(). The optimization method is Newton-Raphson with
step-halving, or the Levenberg-Marquardt method. The latter has
been shown to converge better when there are large offsets, or when
the likelihood surface is moderately ill-conditioned (see
Details). Execution time is fast even for hundreds of thousands
of intercepts. The limiting factor is the number of intercepts
times the number of columns of x.
orm.fit also fits a single-level random-intercept extension of
this model when cluster is specified, for clustered or repeated-
measures data, optionally with a user-specified, continuously- or
discretely-varying weight on the random effect (see mre below and
Details). See Details for a full description of the
random-intercept model, how it is fit, and how sparseness is
preserved even with thousands of intercepts.
orm.fit(
x = NULL,
y,
cluster = NULL,
family = c("logistic", "probit", "loglog", "cloglog", "cauchit"),
offset,
initial,
opt_method = c("NR", "LM"),
maxit = 30L,
eps = 5e-04,
gradtol = 0.001,
abstol = 1e+10,
minstepsize = 0.01,
maxetastep = 15,
tol = .Machine$double.eps,
trace = FALSE,
penalty.matrix = NULL,
weights = NULL,
normwt = FALSE,
scale = FALSE,
mscore = FALSE,
inclpen = TRUE,
y.precision = 7,
compstats = TRUE,
onlydata = FALSE,
sigma.init = 1,
maxit.outer = 100L,
maxit.mode = 30L,
nAGQ = 7L,
nAGQ.grid = c(7L, 11L, 15L, 21L, 31L, 45L, 63L),
nAGQ.tol = 1e-05,
mre = NULL,
...
)a list with the following components, not counting all the
components produced by orm.fit:
call: calling expression
freq: table of frequencies for y in order of increasing y
yunique: vector of sorted unique values of y
stats: vector with the following elements: number of
observations used in the fit, effective sample size ("ESS"; see
Details for a limitation when cluster is given), number of
unique y values, median y from among the observations used in
the fit, maximum absolute value of first derivative of log
likelihood, model likelihood ratio chi-square, d.f., P-value,
score chi-square and its P-value, Spearman's \(\rho\) rank
correlation between linear predictor and y (if there is no
censoring), Somers' \(Dxy\) rank correlation (if there is no
censoring or only right censoring), the Nagelkerke \(R^2\)
index, other \(R^2\) measures, the \(g\)-index, gr (the
\(g\)-index on the ratio scale), and pdm (the mean absolute
difference between 0.5 and the estimated probability that
\(y \geq\) the marginal median). When penalty.matrix is
present, the chi-square, d.f., and P-value are not corrected for
the effective d.f. When cluster is given, the score-test-related
elements are not currently available (see Details) and come
back as NA.
fail: set to TRUE if convergence failed (and maxit>1, or,
when cluster is given, maxit.outer>1)
coefficients: estimated parameters (intercepts then slopes)
family, famfunctions: see orm()
deviance: named vector of -2 log likelihoods, in the order
computed. Element names, in order (present or absent according to
which of censoring/an offset/covariates/clustering actually
apply -- see Details for why, when cluster is given and
p>0, two additional, clustered elements are always computed):
"intercepts" (always present); "intercepts+offset" (only if
an offset variable is given); "intercepts+x" (only if p>0
covariates are present); "intercepts+random effects" (only if
cluster is given -- when p>0 this is an extra, covariate-free
clustered fit computed solely to give the next element a correct
baseline; when p==0 it is instead the final element, there
being no x to add); "intercepts+x+random effects" (the final
element whenever both cluster is given and p>0). The vector
is correspondingly shorter whenever cluster is not given,
ending at "intercepts+x" (or just "intercepts" if p==0).
lpe: vector of per-observation likelihood probability elements.
An observation's contribution to the log likelihood is the log of
lpe. Not currently computed when cluster is given (returned
as NA; see Details).
non.slopes: number of intercepts in model
interceptRef: the index of the middle (median) intercept used
in computing the linear predictor and var
linear.predictors: the linear predictor using the first
intercept. When cluster is given, this is the fixed-effects-only
linear predictor -- random effects are not added in, so that
downstream fit statistics retain their ordinary interpretation.
penalty.matrix: see above
info.matrix: see orm() and infoMxop(). When cluster is
given, this is the corrected (Louis's-identity-adjusted) sparse
information matrix described in Details, with additional
row(s)/column(s) in its b/ab pieces: a single log(sigma)
element when mre is not given, or log(sigma) and sigma2
together when it is. Use infoMxop(fit$info.matrix, i='sigma_parameters') to retrieve the variance of log(sigma)
alone, or the 2x2 covariance matrix of log(sigma) and
sigma2 when both are present.
ncluster: number of distinct clusters. Present only when
cluster is given.
nAGQ: the number of adaptive Gauss-Hermite quadrature nodes
used for the final, converged fit (after any automatic
escalation). Present only when cluster is given.
sigma: the estimated random-intercept standard deviation.
Present only when cluster is given and mre is not; when
mre is given, sigma1 and sigma2 are reported instead (see
below) and sigma is absent.
sigma1, sigma2: the two estimated random-effect scale
parameters described in Details (sigma1 constrained
positive, sigma2 unconstrained). Present only when both
cluster and mre are given.
gamma: vector of estimated (posterior mode) random intercepts,
one per cluster, in the order of the distinct levels of
cluster. Present only when cluster is given.
design matrix with no column for an intercept
response vector, numeric, factor, or character. The
ordering of levels is assumed from factor(y).
optional vector, the same length as y, of cluster
identifiers (e.g. subject ID for repeated measurements). When
given, orm.fit fits an additional normally-distributed random
intercept for each distinct cluster; see Details. Left as
NULL (the default) to fit the ordinary fixed-effects model.
a character value specifying the distribution family,
corresponding to logistic (the default), Gaussian, Cauchy, Gumbel
maximum (exp(-exp(-x)); extreme value type I), and Gumbel
minimum (1-exp(-exp(x))) distributions. These are the
cumulative distribution functions assumed for
Prob[Y >= y | X]. The family argument can be an unquoted or a
quoted string, e.g. family=loglog or family="loglog". To use
a built-in family, the string must be one of the following,
corresponding to the previous list: logistic, probit,
loglog, cloglog, cauchit.
optional numeric vector containing an offset on the logit scale
vector of initial parameter estimates, beginning
with the intercepts. If initial is not specified, the function
computes the overall score \(\chi^2\) test for the global null
hypothesis of no regression. initial is padded to the right
with zeros for the regression coefficients, if needed. When
censoring is present, initial can also be a list with elements
time and surv from the npsurv attribute of the y element
of a previous fit. This is useful when bootstrapping, for
example. When cluster is given, orm.fit's own fixed-effects
initial-value refinement (run once, regardless of cluster) is
reused as the warm start for the random-intercept fit, so
initial need not account for random effects.
set to "LM" to use Levenberg-Marquardt instead
of the default Newton-Raphson; see Details for the asymmetric
damping schedule used.
maximum no. iterations (default=30). For the
random-intercept outer loop this is instead governed by
maxit.outer.
difference in \(-2 \log\) likelihood for declaring
convergence. Default is .0005. This handles the case where the
initial estimates are MLEs, to prevent endless step-halving.
maximum absolute gradient before convergence can be
declared. For the fixed-effects fit, gradtol is automatically
scaled by n / 1000 since the gradient is proportional to the
sample size; for the random-intercept outer loop it is used as an
absolute tolerance instead (see Details).
maximum absolute change in parameter estimates from
one iteration to the next before convergence can be declared; by
default has no effect for the fixed-effects fit. For the
random-intercept outer loop, this is capped at 1e-4 internally
regardless of the value supplied, since the default of 1e10 was
found to make this criterion a permanent no-op for that loop's
only-linear convergence.
used to specify when to abandon step-halving
maximum change, in linear-predictor (eta) units,
that a single Newton-Raphson (opt_method='NR', the default) or
Levenberg-Marquardt (opt_method='LM') update may impose before it
is accepted; a larger implied change forces additional step-halving
('NR') or an outright-rejected, more heavily damped correction
('LM'), the same as the existing intercept-ordering safeguard.
Also governs the analogous joint Newton step of the clustered/
random-intercept fit (used when cluster is given), which shares
the identical risk. Guards against Matrix::solve()'s native
sparse LU factorization crashing on the severely, badly-scaled
Hessian an undamped early step can produce with a very large
number of intercepts (e.g. a nearly continuous y); see
ormfit's internal documentation for the mechanism. The default,
15, is comfortably above the largest step a well-posed fit
legitimately needs and should rarely bind; lower it only if you
continue to see native crashes, and raise it only if you have
reason to believe legitimate large early steps are being
needlessly halved for your data.
singularity criterion. Default is typically 2e-16
set to TRUE to print -2 log likelihood,
step-halving fraction, change in -2 log likelihood, maximum
absolute value of first derivative, and max absolute change in
parameter estimates at each iteration. For the random-intercept
outer loop, set to a value greater than 1 to also print
per-outer-iteration mode-finding and AGQ diagnostics.
a self-contained ready-to-use penalty matrix
-- see lrm
a vector (same length as y) of possibly fractional
case weights
set to TRUE to scale weights so they sum to n,
the length of y; useful for sample surveys as opposed to the
default of frequency weighting
set to TRUE to subtract column means and divide by
column standard deviations of x before fitting, and to
back-solve for the un-normalized covariance matrix and regression
coefficients. This can sometimes make the model converge for very
large sample sizes where for example spline or polynomial
component variables create scaling problems leading to loss of
precision when accumulating sums of squares and crossproducts.
set to TRUE to compute the sparse score matrix and
store its elements as a list mscore. Not currently computed
when cluster is given.
set to FALSE to not include the penalty matrix in
the Hessian when the Hessian is being computed on transformed
x, vs. adding the penalty after back-transforming. This should
not matter.
when y is numeric, values may need to be
rounded to avoid unpredictable behavior with unique() with
floating-point numbers. Default is to 7 decimal places.
set to FALSE to prevent the calculation of the
vector of model statistics
set to TRUE to return the data used in model
fitting as a list, without fitting the model
starting value for the random-intercept standard
deviation sigma, used only when cluster is given
maximum number of outer-loop iterations for the
random-intercept fit (mode-finding / AGQ-reweighted Newton step /
sigma profile search), used only when cluster is given
maximum number of Newton-Raphson iterations for
the per-cluster random-intercept mode-finding inner loop, used
only when cluster is given
starting number of adaptive Gauss-Hermite quadrature
nodes for the random-intercept fit, used only when cluster is
given; escalated automatically as needed (see Details)
increasing sequence of candidate quadrature-node
counts nAGQ is escalated through, used only when cluster is
given
relative log-likelihood change, comparing the
current nAGQ to the next candidate in nAGQ.grid, below which
nAGQ escalation stops; used only when cluster is given
optional fixed, numeric vector the same length as y,
used only when cluster is also given: a user-computed
"multiplier for random effects" governing how the random
intercept's own contribution is weighted at each observation, via
two estimated scale parameters sigma1 and sigma2 in place of
the ordinary single sigma (see Details). Left as NULL
(the default) to fit the ordinary, single-sigma random-
intercept model. For a Markov-1 model on repeated ordinal
responses, a natural choice is mre equal to 0 at each
subject's first follow-up and 1 at every later one (see
Examples), but any fixed function of the data, continuous or
discrete, is valid. mre must vary within at least some
clusters; a completely constant mre (whether 0, 1, or any
other single value across the whole dataset) is rejected with an
informative error, since it would leave sigma1 and sigma2
unidentifiable.
ignored
Frank Harrell
Department of Biostatistics, Vanderbilt University
[email protected]
When cluster is given, orm.fit fits k ordinal intercepts and
p regression coefficients as usual, plus a normally-distributed
random intercept gamma_j ~ N(0, sigma^2) for each of the distinct
clusters (levels of cluster). Fitting integrates out the random
effects to obtain the marginal likelihood, using adaptive
Gauss-Hermite quadrature (AGQ) rather than a Laplace approximation.
This choice was empirically driven: extensive validation against
ordinal::clmm2 showed that Laplace approximation (equivalent to
nAGQ=1) produces substantial, systematic bias in both the fixed
effects and, especially, the variance component -- worsening as
cluster informativeness decreases (small clusters, large sigma),
which is the regime most common in practice. nAGQ is not held
fixed either: it is escalated automatically, refitting (warm-
started) at each successively larger value in nAGQ.grid until the
log-likelihood stabilizes within nAGQ.tol, since no single fixed
value of nAGQ was found to be adequate across the full range of
sigma a user might encounter.
mre
When a fixed, user-specified per-observation vector mre is also
given, the random intercept's own contribution to the linear
predictor is no longer a constant gamma_j -- it is scaled by a
quantity that can vary, continuously if desired, across an
individual cluster's own observations:
$$\eta_{ij} = \ldots + [\sigma_1(1-\mbox{mre}_{ij}) +
\sigma_2\,\mbox{mre}_{ij}]\, v_i$$
where v_i ~ N(0,1) is a standardized, per-cluster random effect
and sigma1, sigma2 are two estimated scale parameters --
sigma1 constrained positive via a log link (to pin down v_i's
otherwise arbitrary sign, since flipping the sign of both sigma2
and every v_i simultaneously leaves the likelihood unchanged
unless something anchors it), sigma2 left unconstrained, so it
may legitimately come out negative. mre is entirely user-computed
and fixed -- it is data, exactly like a column of x, never
estimated -- while sigma1 and sigma2 are estimated jointly with
alpha and beta in the same Newton step. mre=NULL (the
default) recovers the ordinary, single-sigma model exactly, with
no loss of speed or accuracy.
This additive form -- sigma1 and sigma2 each multiplying a
separately-known quantity, (1-mre)v_i and mre \, v_i
respectively -- was chosen over an earlier, multiplicative design
(a single sigma combined with a separately estimated multiplier)
specifically because it keeps sigma1 and sigma2 close to
orthogonal. Extensive testing found the multiplicative form's extra
parameter became structurally unidentifiable whenever the base
sigma approached its lower boundary of zero, causing the fit to
converge to an arbitrary, meaningless value for it; sigma1 and
sigma2 remain well-behaved instead (their estimated correlation
stayed well under 0.2 in magnitude across every scenario tested,
including deliberately small samples and scale parameters of
opposite sign), since each is an absolute scale directly identified
by its own share of the data rather than a ratio relative to the
other.
mre must vary within at least some clusters for sigma1 and
sigma2 to be separately identifiable; orm.fit checks this
directly and stops with an informative error if mre has no
variation anywhere in the data (e.g. is constant at 0, at 1, or
at any other single value for every observation). A cluster
contributing only a single observation is harmless and needs no
special handling -- it simply contributes nothing toward separating
sigma1 from sigma2, the same way a covariate with no variation
in one subgroup does not hurt an ordinary regression, as long as
some clusters elsewhere in the data do have within-cluster mre
variation.
A natural and recommended choice for a Markov-1 model on repeated
ordinal responses (a proportional-odds model whose linear predictor
includes a term for the previous response) is a piecewise-constant
mre: 0 at each subject's first follow-up time and 1 at every
subsequent one (see Examples). The motivation is that a plain,
constantly-weighted random effect double-counts its own influence
once a lag term is present -- the previous response already carries
the random effect's influence forward, so adding it again at full
strength at every subsequent visit can induce an unrealistically
large or ever-growing correlation across time. Letting sigma2
differ from sigma1 -- including, if the data call for it, being
of the opposite sign -- lets the random effect's influence after
the first visit either add to or subtract from the correlation the
lag term alone would induce, as the data determine. mre is not
restricted to this 0/1 step function, however; any fixed,
user-computed function of the data (continuously varying with
follow-up time, for example) is equally valid.
When cluster is given together with p>0 covariates, orm.fit
fits an additional, otherwise-unnecessary clustered model:
intercepts plus the same random-effects structure (including mre,
if given), but without x. This provides a clustered null deviance
that is directly comparable to the clustered full model's deviance,
so that the reported model likelihood ratio chi-square is a genuine
test of the covariates given the random-effects structure already
in the model -- not a test confounded with whether the random
effect itself improves the fit. See the deviance element of
Returns for the full, named composition of the returned
deviance vector.
Fitting uses an outer, block-alternating loop, each iteration of which does the following:
Per-cluster mode-finding (the inner loop): a scalar
Newton-Raphson iteration finds each cluster's own posterior mode
and curvature for its random intercept, at the current
(alpha, beta, sigma). This is vectorized across all clusters
simultaneously -- one lean Fortran call per Newton step, not one
call per cluster -- so its cost is O(n) regardless of the
number of clusters.
AGQ-reweighted Newton step for (alpha, beta): at each of
the current nAGQ quadrature nodes (centered and scaled per
cluster using the mode and curvature just found), the same,
completely unmodified orm.fit likelihood/gradient/Hessian
machinery is called once, with that node's random-intercept
value folded into the offset and its posterior node weight
folded into the case weights. The nAGQ results are summed to
form one Newton step. Because this reuses orm.fit's own
tri-band (or, under interval censoring, general sparse triplet)
intercept-Hessian representation and its O(k) solve unchanged,
the random-intercept extension inherits the same scalability to
very large numbers of intercepts that orm.fit has without
random effects -- it simply pays this cost nAGQ times per
outer iteration instead of once.
Fitting the random-effect scale parameter(s). When mre is
not given, this is a one-dimensional profile search for sigma
via optimize(), using only the lean, O(n) per-observation
routine from step 1 (no gradient or Hessian needed) -- cheap
relative to step 2. This separate-search design was tried for
the sigma1/sigma2 case too, but found to converge to a
persistent wrong value whenever the true sigma2 had the
opposite sign from its (necessarily sign-unknown-in-advance)
starting value, regardless of sample size -- confirmed to be a
genuine optimization-path problem, not a sign of poor
identifiability, since a direct grid-scan of the likelihood
always found a well-behaved optimum near the true values. When
mre is given, sigma1 and sigma2 are therefore instead
estimated jointly with alpha and beta in the very same
Newton step as step 2, sigma2 started at exactly 0 (not an
arbitrary nonzero value) so the step itself discovers the
correct sign directly from the data rather than being biased
toward whichever sign it happened to start on.
Because this scheme converges only linearly near the optimum
(unlike orm.fit's own Newton-Raphson, which converges
quadratically), its convergence criteria differ correspondingly:
the gradient tolerance is used as an absolute, not an n-scaled,
criterion (an n-scaled tolerance is nearly a no-op at realistic
sample sizes for a scheme with only linear convergence), and a
separate allowance accepts convergence once sigma's absolute
change falls below a small floor even if its relative change has
not, since relative changes become numerically meaningless once
sigma is already effectively at its lower boundary of zero.
The outer loop's own convergence is accelerated using SQUAREM (Varadhan and Roland, 2008), a generic extrapolation method for fixed-point iterations that share this scheme's linear/geometric convergence pattern. Every SQUAREM cycle runs two ordinary outer iterations, uses them to estimate the local geometric contraction rate, and extrapolates a candidate point that would take many more ordinary iterations to reach directly. That candidate is never accepted on faith: it is itself run through one more ordinary outer iteration ("stabilization"), and is kept only if that step's own parameter change is no larger than the plain, unaccelerated iteration's would have been -- so acceleration can only match or improve on plain iteration, never meaningfully underperform it. The extrapolation's step length is additionally bounded by an adaptive cap that starts conservative and grows after each accepted extrapolation, shrinking back after a rejected one -- an unbounded step length is only trustworthy once the local trajectory has settled into a clean geometric pattern, and can otherwise produce wild, wastefully expensive trial points. In testing this combination reduced wall clock fitting time by roughly 15-20% relative to unbounded extrapolation, and reduced the number of outer iterations needed by roughly 3-to-7-fold relative to no acceleration at all, on both a very small (20-observation, 10-cluster) and a more typical-scale (400-observation, 200-cluster) dataset, while landing at coefficient estimates matching plain (unaccelerated) iteration to well within the convergence tolerance itself.
Every piece of the random-intercept computation was designed so
that its cost scales with n, the number of clusters, or cluster
size -- never with k^2 or k^3:
The per-cluster mode-finding and AGQ node evaluation touch k
only as a lookup into alpha, never as a matrix dimension.
The (alpha, beta) Newton step reuses orm.fit's existing
tri-band (no interval censoring) or general sparse triplet
(interval censoring present) representation of the intercept
Hessian and its O(k) tridiagonal solve, completely unchanged.
The final, corrected information matrix's missing-information
correction (via Louis's (1982) identity, which corrects for the
variance a fixed-posterior-weight Hessian omits) is accumulated
as a sparse triplet list: each cluster contributes only to the
alpha indices that cluster's own observations actually touch --
at most about twice its own size, however large k is. The
total number of nonzero entries this produces is bounded by
roughly the number of clusters times the square of a typical
cluster's size, independent of k, for realistic (bounded
cluster size) designs -- confirmed empirically (e.g. 588 nonzero
entries out of 11,175 possible for a k=149, 75-cluster
example). This lets orm.fit fit random-intercept models with
thousands of intercepts in about a second, versus tens of seconds
to minutes for ordinal::clmm2 on the same data.
log(sigma)'s standard error (or the covariance of log(sigma)
and sigma2, when mre is given) is obtained by extending only
the non-intercept (b/ab) pieces of info.matrix by one or two
rows and columns -- the intercept-intercept block is never
touched by this extension, so infoMxop() and other code that
consumes info.matrix require no special handling for the
random-effects case beyond what interval censoring already
required. infoMxop() locates these elements by searching
info.matrix$xname for "log(sigma)" and "sigma2" wherever
they occur, rather than assuming a fixed position, so their
presence and order elsewhere in xname does not matter.
opt_method='LM' uses an asymmetric Levenberg-Marquardt damping
schedule (lambda divided by 1.5 on an accepted step, multiplied
by 10 on a rejected one) rather than a symmetric /10/*10 rule.
The symmetric rule can enter a permanent limit cycle on
moderately ill-conditioned (but not literally singular)
likelihood surfaces -- alternating forever between a damping
level that barely succeeds and one that barely fails, with the
true optimum never actually reached. The asymmetric rule was
confirmed, by direct inspection of lambda's trajectory on such
a case, to break this cycle, at no cost in iterations for
well-conditioned problems.
Fitted standard errors (via infoMxop() operating on
info.matrix) always come from this corrected, not naive,
information matrix when random effects are present -- ignoring
the missing-information correction was found in validation
studies to understate intercept standard errors by as much as
40%.
Requesting i='sigma_parameters' from infoMxop() can encounter
a genuinely singular information matrix when a fit has converged
with sigma (or sigma1) at or very near its lower boundary of
zero -- at that boundary its information content correctly
collapses to zero (a known, non-regular feature of variance-
component estimation, not a numerical defect), and infoMxop
returns NA (a scalar, or an all-NA 2x2 matrix when sigma2
is also present) with an explanatory warning rather than raising
a raw singularity error in that case.
The effective-sample-size adjustment for censored observations
(contributing to the "ESS" element of stats) is not yet
available when cluster is specified; stats["ESS"] reflects
only the uncensored-observation count in that case, and a warning
is issued.
orm(), lrm(), glm(), gIndex(), SparseM::solve(),
recode2integer(), infoMxop()
if (FALSE) {
# Fit an additive logistic model containing numeric predictors age,
# blood.pressure, and sex, assumed to be already properly coded and
# transformed
#
# fit <- orm.fit(cbind(age, blood.pressure, sex), death)
# Fit a random-intercept version for repeated measurements on
# subject id
# fit <- orm.fit(cbind(age, blood.pressure, sex), death, cluster=id)
# A Markov-1 model on repeated ordinal responses, where the linear
# predictor already includes a term for the previous response (e.g.
# via an added column of x, such as a spline in the lagged y). The
# recommended mre for this situation takes on only the values 0 and
# 1: 0 at each subject's FIRST follow-up time, 1 at every later one.
# This lets the fit determine, via the separately-estimated sigma2,
# how much the random effect should add to or subtract from the
# correlation the lag term alone induces, rather than assuming (as
# a plain, constantly-weighted random effect implicitly would) that
# its influence should be injected again at full, undiminished
# strength at every subsequent visit:
# mre <- ave(time, id, FUN = function(t) as.numeric(t > min(t)))
# fit <- orm.fit(cbind(prev.y, age, blood.pressure, sex), y,
# cluster=id, mre=mre)
}
Run the code above in your browser using DataLab