Processes four types of information matrices: ones produced by the SparseM package for the orm or lrm functions in rms version 6.9-0 and earlier, by the Matrix package for version 7.0-0 of rms using a tri-band diagonal matrix for the intercepts, using Matrix for general sparse information matrices for intercepts (when any interval-censored observations exist or random effects are included), or plain matrices. For Matrix, the input information matrix is a list with three elements: a containing in two columns the diagonal and superdiagonal for intercepts (when there is no interval censoring) or a list with three elements row, col, a (when there is interval censoring), b, a square matrix for the covariates, and ab for intercepts x covariates. If nothing else is specified, the assembled information matrix is returned for Matrix, or the original info otherwise. If p=TRUE, the number of parameters in the model (number of rows and columns in the whole information matrix minus any log(sigma)/sigma2 elements) is returned. If i is given, the i elements of the inverse of info are returned, using efficient calculation to avoid inverting the whole matrix. Otherwise if invert=TRUE or B is given without i, the efficiently (if Matrix or SparseM) inverted matrix less any log(sigma)/sigma2 elements is returned, or the matrix multiplication of the inverse and B. If both i and B are given, what is returned is the i portion of the inverse of the information matrix, matrix multiplied by B. This is done inside solve().
infoMxop(
info,
i,
invert = !missing(i) || !missing(B),
B,
np = FALSE,
tol = .Machine$double.eps,
abort = TRUE
)a single integer or a matrix
an information matrix object
integer vector specifying elements returned from the inverse. You can also specify i='x' to return non-intercepts, i='i' to return intercepts, or i="sigma_parameters" to return the estimated variance of log(sigma) (a scalar), or the 2x2 covariance matrix of log(sigma) and sigma2 when both are present (a two-scale random-effect model).
set to TRUE to invert info (implied when i or B is given)
multiplier matrix. If random effects were included, B should ignore them.
set to TRUE to just fetch the total number of parameters (intercepts + betas, any log(sigma)/sigma2 not counting)
tolerance for matrix inversion singularity
set to FALSE to run the solve calculation through try() without aborting; the user will detect that the operation did not success by examinine inherits(result, 'try-error') for being TRUE. This is ignored when i='sigma_parameters': if the underlying solve() fails there (e.g. because the model converged with a scale parameter at or near its lower boundary of zero, where it is not identifiable), a warning is issued and NA is returned regardless of abort, since a plain NA (scalar or all-NA 2x2 matrix, as appropriate) is more informative than a try-error object for this specific, well-understood failure mode. For every other request, a singularity caused specifically by one of these scale parameters being at this same boundary is instead handled by an automatic, targeted ridge retry (see Description) with a warning on success; abort only comes into play if that retry itself fails, in which case the original error is handled exactly as documented above.
Frank Harrell
When only variance-covariance matrix elements corresponding to the non-intercepts are desired, specify
i='x' or i=(k + 1) : nv where nv is the number of intercepts and slopes combined. infoMxop computes the needed covariance matrix very quickly in this case. To retrieve only the variance/covariance of the random-effect scale parameter(s) for a model with random effects, specify i='sigma_parameters'.
When inverting info, if info has a 'scale' attribute with elements mean and sd, the scaling is reversed after inverting info.
When the number of intercepts is needed to be known and the info object is not a 3-element list, info must have an intercepts attribute to define the number. This is used for example when transx=TRUE is specified to lrm or lrm.fit.
When the model contains random effects, the xname component of the info list must contain an element named log(sigma) (the log of the random effect's scale, or of its first scale parameter in a two-scale model) and, for a two-scale (sigma1/sigma2) random-effect model, an additional element named sigma2. infoMxop locates these BY NAME within xname, wherever they occur -- it does not assume they are the last element(s) of xname or in any particular order relative to each other.
A random-intercept fit's information matrix can become genuinely near-singular in one or two directions -- log(sigma), and sigma2 when present -- whenever the corresponding scale parameter is at or very near a boundary where it is not identifiable; this is a real, expected non-regular feature of variance-component estimation, not a numerical defect. Since every general solve (missing(i), i='i', or any i mixing intercepts and slopes) factors the whole information matrix regardless of which elements of its inverse are ultimately kept, this near-singularity can surface even when the caller never asked for sigma_parameters at all. When it does, infoMxop automatically retries with a small, escalating ridge applied only to the affected diagonal entry/entries -- not the whole matrix -- since these are empirically essentially uncorrelated with everything else at such a fit, so this leaves every other requested quantity negligibly perturbed. A warning is issued when this happens. i='x' is never affected by this, since it already avoids factoring the full matrix (see below); i='sigma_parameters' is also never affected, since that variance/covariance genuinely isn't well defined there and is instead returned as NA (a scalar or an all-NA 2x2 matrix, matching whether one or two scale parameters are present) with its own explanatory warning.
if (FALSE) {
f <- orm(y ~ x)
infoMxop(f$info.matrix) # assembles 3 pieces
infoMxop(v, i=c(2,4)) # returns a submatrix of v inverse
infoMxop(f$info.matrix, i='x') # sub-covariance matrix for just the betas
f <- orm(y ~ x + cluster(id))
infoMxop(f$info.matrix, i='sigma_parameters')
# a single number = variance of log(sigma), sigma=sqrt of random effects variance,
# OR (for a two-scale sigma1/sigma2 random-effect model) a 2x2 covariance matrix
# for log(sigma) and sigma2
}
Run the code above in your browser using DataLab