Skip to contents

Compare Standard Errors from orm and vglm

Source: R/robcov_vglm.R
compare_se_orm_vglm.Rd

Utility function to compare standard errors between rms::orm() and VGAM::vglm() fits, including robust (sandwich) versions. This is useful for validating that both approaches give similar results.

Usage

compare_se_orm_vglm(orm_fit, vglm_fit, cluster = NULL)

Arguments

orm_fit

A fitted orm object from the rms package.

vglm_fit

A fitted vglm object from the VGAM package. Should be fit with cumulative(parallel = TRUE, reverse = TRUE) to match orm's parameterization for ordinal outcomes, or binomialff for binary outcomes.

cluster

Optional cluster variable for robust standard errors. Must have the same length as the number of observations in both models.

Value

A data frame comparing standard errors with columns:

parameter

Parameter name (intercepts numbered, then regression coefficients)

se_orm_model

Model-based SE from orm

se_vglm_model

Model-based SE from vglm

se_orm_robust

Robust SE from orm

se_vglm_robust

Robust SE from vglm

ratio_model

Ratio of vglm to orm model-based SEs

ratio_robust

Ratio of vglm to orm robust SEs

Details

This function aligns parameters between the two model types:

  • orm uses names like "y>=2", "y>=3", etc. for intercepts (or "Intercept" for binary outcomes)

  • vglm uses names like "(Intercept):1", "(Intercept):2", etc.

The comparison focuses on:

  1. Model-based standard errors (from the inverse Hessian)

  2. Robust (sandwich) standard errors

Ratios close to 1.0 indicate good agreement. For equivalent cumulative-logit likelihoods, the default observed-score robust standard errors should agree closely after parameter alignment. Native model-based standard errors can still differ because VGAM reports final IRLS working-information covariance.

Examples

if (FALSE) { # \dontrun{
library(rms)
library(VGAM)

# Fit models
dd <- datadist(mydata)
options(datadist = "dd")
m_orm <- orm(y ~ x1 + x2, data = mydata, x = TRUE, y = TRUE)
m_vglm <- vglm(y ~ x1 + x2,
               family = cumulative(parallel = TRUE, reverse = TRUE),
               data = mydata)

# Compare standard errors
comparison <- compare_se_orm_vglm(m_orm, m_vglm)
print(comparison)

# With clustering
comparison_cl <- compare_se_orm_vglm(m_orm, m_vglm, cluster = mydata$id)
} # }