Joint Covariance Matrix for gamlss2 Models

Description

Computes coefficient covariance from the joint penalized information of a fitted gamlss2 model.

Usage

## S3 method for class 'gamlss2'
vcov(object,
  type = c("vcov", "cor", "se", "coef"), full = FALSE,
  method = c("joint", "working", "numeric"), unconditional = FALSE, ...)

Arguments

object A fitted gamlss2 model.
type The covariance matrix, correlation matrix, standard errors, or corresponding coefficient vector to return.
full Logical. If FALSE, return ordinary linear coefficients. Their covariance is still extracted from the complete joint inverse. If TRUE, include supported smooth coefficients.
method Information calculation. “joint” uses the raw joint observed information, “working” uses the final family working curvature with zero cross-parameter blocks, and “numeric” uses a slow coefficient-dimensional numerical Hessian for validation.
unconditional Logical. If TRUE, a correction for uncertainty in selected smoothness parameters is added. For JR fits this uses the fitted joint REML curvature. Other fits use a post-fit local approximation. This requires method = “joint”.
… Control arguments passed to the numerical Hessian when method = “numeric”.

Details

The default covariance is the inverse of \(H + P\), where \(H\) is the joint likelihood information for every supported distributional predictor and \(P\) is the complete penalty at the fitted smoothing and regularization parameters. Thus it retains covariance among linear and smooth terms and among distributional parameters. Prior observation weights and offsets in the final linear predictors are included.

Joint information is computed by differentiating the linked family scores. Mixed derivatives are checked in both orientations. If finite differences disagree, additional step sizes and Richardson extrapolation are used to check numerical convergence; persistent inconsistencies raise an error.

For monotonic ms terms, inversion is restricted to coefficient directions that preserve the fitted active constraints and term centering. This local approximation conditions on the active constraint set and does not include uncertainty about which constraints are active. Unconditional and sandwich corrections are currently unavailable for these terms.

The default covariance conditions on fitted smoothing parameters. Use unconditional = TRUE to include their approximate estimation uncertainty. Fixed smoothing parameters receive no correction, and model-selection uncertainty is excluded. Adaptive smooth penalties, selected ridge penalties, and non-estimable coefficients are unsupported.

Special terms that have no verifiable coefficient-linear design are rejected explicitly. Objects fitted with light = TRUE do not retain the final predictors and must be refitted for joint covariance. Aliased or numerically unidentified coefficient variances are returned as NA.

The ordinary summary() and effect-result paths use the covariance stored with each fitted model term by default, which is substantially faster for large models. Use summary(object, method = “joint”) or results(object, method = “joint”) to request the complete joint information. Use summary(object, unconditional = TRUE) for the smoothness correction. This option selects the joint method unless another method is specified explicitly. Simulation-based predict() calls always use the complete joint covariance. No consumer silently substitutes method = “working” when observed joint information is indefinite.

Value

A matrix for type = “vcov” or “cor”; a named numeric vector for type = “se” or “coef”.

See Also

gamlss2, predict.gamlss2

Examples

library("gamlss2")

data("cars", package = "datasets")
m <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)

## linear-coefficient covariance,
## marginalized over smooth coefficients.
vcov(m)

## complete linear and smooth covariance
V <- vcov(m, full = TRUE)
head(diag(V))

## include estimated smoothness uncertainty
Vc <- vcov(m, full = TRUE, unconditional = TRUE)