Robust (Sandwich) Covariance Matrix Estimation for vglm Objects

Computes the Huber-White sandwich covariance matrix estimator for vglm objects, optionally with cluster-robust standard errors. This provides functionality for vglm models using the same HC1 and finite-cluster scalar corrections as the package-owned ORM sandwich.

Usage

robcov_vglm(
  fit,
  cluster = NULL,
  adjust = NULL,
  bread = c("observed", "vglm"),
  type = c("HC0", "HC1"),
  cadjust = NULL
)

Arguments

fit: A fitted vglm object from the VGAM package.
cluster: Optional vector of cluster identifiers. If provided, computes cluster-robust standard errors. Must have the same length as the fitted observations, or the same length as the original pre-NA data when the fitted model records omitted rows in fit@na.action.
adjust: Deprecated compatibility alias for cadjust. A non-missing logical value is accepted so existing calls continue to work. Do not supply both adjust and cadjust.
bread: Character string selecting the sandwich bread. "observed" (the default) numerically differentiates the summed analytic scores. "vglm" uses the inverse final VGAM IRLS working information returned by stats::vcov(fit).
type: Character string selecting an HC correction. "HC0" applies no observation degrees-of-freedom correction. "HC1" multiplies the meat by $(m - 1) / (m - p)$, where $m$ is the number of fitting rows and $p$ is the number of coefficients.
cadjust: Logical. Apply the cluster correction $n / (n - 1)$, where $n$ is the number of clusters? The default is TRUE when cluster is supplied and FALSE otherwise. When explicitly TRUE without cluster, each fitted row is treated as its own cluster for this scalar correction.

Value

A list with class "robcov_vglm" containing all original model information plus robust covariance estimates:

Model information (copied from original vglm fit):

coefficients: Model coefficients
fitted.values: Fitted values from the model
residuals: Model residuals
family: The VGAM family object
predictors: Linear predictors (eta)
prior.weights: Prior weights
df.residual: Residual degrees of freedom
df.total: Total degrees of freedom
rank: Model rank
criterion: Convergence criterion values
constraints: Constraint matrices
control: Control parameters
extra: Extra information from family
misc: Miscellaneous model info
model: Model frame (if available)
x: Design matrix (if available)
y: Response (if available)
terms: Model terms
assign: Assignment vector
xlevels: Factor levels
na.action: NA handling
call: Original model call
original_call: Call to robcov_vglm

Robust covariance estimates:

var: The robust variance-covariance matrix
se: Robust standard errors (sqrt of diagonal)
z: Z-statistics (coefficients / robust SE)
pvalues: Two-sided p-values based on z-statistics
original_vcov: The original model-based variance-covariance matrix
original_se: Original model-based standard errors

Sandwich components:

scores: Matrix of observation-level score contributions ($m \times p$)
influence_beta_obs: Matrix of observation-level influence contributions for the coefficient vector $\beta$ ($m \times p$), computed as the scores times the selected bread
bread: The selected covariance-scale bread
bread_type: The requested bread type
meat: The crossproduct of raw observation scores, or cluster-summed scores when cluster is supplied

Clustering information:

cluster: The cluster variable (if provided)
n: Number of observations
n_clusters: Number of clusters (if clustered)
type: The HC correction type
cadjust: Whether the cluster correction was applied
adjustment_factor: The combined HC and cluster scalar correction
adjust: Compatibility alias for cadjust

Details

Let $\beta$ be the vector of all $p$ model coefficients and $\hat\beta$ its fitted value. Before optional scalar corrections, the estimated coefficient covariance has the form: $$V_\beta = B^\top M B$$ where $B$ is the selected inverse total information (the sandwich "bread") and $M$ is the score crossproduct (the "meat"). Both are $p \times p$ matrices, and $\top$ denotes transpose. The selected bread is symmetric, so this also equals $B M B$. By default, $B$ is the inverse observed information obtained from the numerical Jacobian of the summed analytic scores. This keeps the meat and bread tied to the same estimating equations and provides misspecification- robust inference. Set bread = "vglm" to use VGAM's final fitted working/Fisher information instead.

Numerical differentiation operates on VGAM's effective VLM coefficients, so proportional-odds, fully nonparallel, and custom partial-proportional-odds constraint matrices use the same implementation. Perturbation steps are reduced adaptively when a candidate produces invalid fitted values.

For unclustered data, the meat matrix is: $$M = \sum_{r=1}^{m} \psi_r \psi_r^\top$$ where $m$ is the number of fitting rows and $\psi_r = \partial\ell_r / \partial\beta$ is the length-$p$ column score vector for row $r$, evaluated at $\hat\beta$. Here $\ell_r$ is that row's log likelihood.

For clustered data, scores are summed within clusters before computing the meat matrix: $$s_i = \sum_{r \in \mathcal{I}_i} \psi_r, \qquad M = \sum_{i=1}^{n} s_i s_i^\top$$ where $i = 1,\ldots,n$ indexes patients (or other independent clusters) and $\mathcal{I}_i$ contains the fitting-row indices for cluster $i$.

The type and cadjust corrections are separate, matching the controls used by sandwich::meatCL(). For clustered data, cadjust = TRUE multiplies the meat by $n / (n - 1)$. type = "HC1" independently multiplies it by $(m - 1) / (m - p)$. Set cadjust = FALSE and type = "HC0" for an unadjusted sandwich.

Z-statistics and p-values: The returned object includes z-statistics computed as coefficients divided by robust standard errors, and two-sided p-values from the standard normal distribution. This approximation is asymptotic in the number of independent observations or clusters. With few clusters, p-values may be anti-conservative.

Weights and aggregated responses: Scores are computed at the fitted-row level. With prior or frequency weights, or with aggregated ordinal count responses, this treats each fitted row as the independent score contribution. That can differ from a sandwich estimator computed after expanding data to one row per independent individual.

Backend comparisons: For equivalent models, results should be compared after accounting for parameterization, cutpoint direction, response coding, weights, bread choice, and missing-data handling.

Examples