Skip to contents

Get variance-covariance matrix from models

Source: R/get_varcov.R
get_varcov.Rd

Returns the variance-covariance, as retrieved by stats::vcov(), but works for more model objects that probably don't provide a vcov()-method.

Usage

get_varcov(x, ...)

# Default S3 method
get_varcov(x, verbose = TRUE, vcov = NULL, vcov_args = NULL, ...)

# S3 method for class 'glmgee'
get_varcov(x, verbose = TRUE, vcov = "robust", ...)

# S3 method for class 'hurdle'
get_varcov(
  x,
  component = "conditional",
  vcov = NULL,
  vcov_args = NULL,
  verbose = TRUE,
  ...
)

# S3 method for class 'aov'
get_varcov(x, complete = FALSE, verbose = TRUE, ...)

# S3 method for class 'mixor'
get_varcov(x, effects = "all", verbose = TRUE, ...)

Arguments

x

A model.

...

Currently not used.

verbose

Toggle warnings.

vcov

Variance-covariance matrix used to compute uncertainty estimates (e.g., for robust standard errors). This argument accepts a covariance matrix, a function which returns a covariance matrix, or a string which identifies the function to be used to compute the covariance matrix.

  • A covariance matrix

  • A function which returns a covariance matrix (e.g., stats::vcov())

  • A string which indicates the kind of uncertainty estimates to return.

    • Heteroskedasticity-consistent: "HC", "HC0", "HC1", "HC2", "HC3", "HC4", "HC4m", "HC5". See ?sandwich::vcovHC

    • Cluster-robust: "CR", "CR0", "CR1", "CR1p", "CR1S", "CR2", "CR3". See ?clubSandwich::vcovCR

    • Bootstrap: "BS", "xy", "residual", "wild", "mammen", "fractional", "jackknife", "norm", "webb". See ?sandwich::vcovBS

    • Other sandwich package functions: "HAC", "PC", "CL", "OPG", "PL".

    • Kenward-Roger approximation: "kenward-roger". See ?pbkrtest::vcovAdj.

    • Finite Population Correction: "fpc" applies the finite population correction. Requires the population_size to be specified in vcov_args. For mixed models, FPC is based on Lai et al. 2018. When vcov = "fpc", at least one of population_size (size of the finite population, must be larger than the number of observations in the model) or cluster_size (the finite size of cluster groups in the population, must be larger than the number of groups of the random effects) in the vcov_args argument. You can additionally apply the Kenward-Roger approximation with the kr argument, e.g. vcov_args = list(cluster_size = 15, kr = TRUE).

Exceptions are following models:

  • Model of class glmgee, which have pre-defined options for the variance-covariance matrix calculation. These are "robust", "df-adjusted", "model", "bias-corrected", and "jackknife". See ?glmtoolbox::vcov.glmgee for details.

  • Model of class glmmTMB currently only support the "HC0" option.

vcov_args

List of arguments to be passed to the function identified by the vcov argument. This function is typically supplied by the sandwich or clubSandwich packages. Please refer to their documentation (e.g., ?sandwich::vcovHAC) to see the list of available arguments. If no estimation type (argument type) is given, the default type for "HC" equals the default from the sandwich package; for type "CR", the default is set to "CR3". For vcov = "fpc", vcov_args must specify either population_size or cluster_size, depending on the model.

component

Should the complete variance-covariance matrix of the model be returned, or only for specific model components only (like count or zero-inflated model parts)? Applies to models with zero-inflated component, or models with precision (e.g. betareg) component. component may be one of "conditional", "zi", "zero-inflated", "dispersion", "precision", or "all". May be abbreviated. Note that the conditional component also refers to the count or mean component - names may differ, depending on the modeling package. See section Model components for details. For models of class glmmTMB, the component argument can also be "full", to return the full variance-covariance matrix (including random effects, called theta).

complete

Logical, if TRUE, for aov, returns the full variance-covariance matrix.

effects

Should the complete variance-covariance matrix of the model be returned, or only for specific model parameters only? Currently only applies to models of class mixor and MixMod.

Value

The variance-covariance matrix, as matrix-object.

Note

get_varcov() tries to return the nearest positive definite matrix in case of negative eigenvalues of the variance-covariance matrix. This ensures that it is still possible, for instance, to calculate standard errors of model parameters. A message is shown when the matrix is negative definite and a corrected matrix is returned.

Model components

Possible values for the component argument depend on the model class. Following are valid options:

  • "all": returns all model components, applies to all models, but will only have an effect for models with more than just the conditional model component.

  • "conditional": only returns the conditional component, i.e. "fixed effects" terms from the model. Will only have an effect for models with more than just the conditional model component.

  • "smooth_terms": returns smooth terms, only applies to GAMs (or similar models that may contain smooth terms).

  • "zero_inflated" (or "zi"): returns the zero-inflation component.

  • "dispersion": returns the dispersion model component. This is common for models with zero-inflation or that can model the dispersion parameter.

  • "instruments": for instrumental-variable or some fixed effects regression, returns the instruments.

  • "nonlinear": for non-linear models (like models of class nlmerMod or nls), returns staring estimates for the nonlinear parameters.

  • "correlation": for models with correlation-component, like gls, the variables used to describe the correlation structure are returned.

  • "location": returns location parameters such as conditional, zero_inflated, smooth_terms, or instruments (everything that are fixed or random effects - depending on the effects argument - but no auxiliary parameters).

  • "distributional" (or "auxiliary"): components like sigma, dispersion, beta or precision (and other auxiliary parameters) are returned.

Special models

Some model classes also allow rather uncommon options. These are:

  • mhurdle: "infrequent_purchase", "ip", and "auxiliary"

  • BGGM: "correlation" and "intercept"

  • BFBayesFactor, glmx: "extra"

  • averaging:"conditional" and "full"

  • mjoint: "survival"

  • mfx: "precision", "marginal"

  • betareg, DirichletRegModel: "precision"

  • mvord: "thresholds" and "correlation"

  • clm2: "scale"

  • selection: "selection", "outcome", and "auxiliary"

  • lcmm: "membership", "longitudinal", "beta", "splines", and "linear"

For models of class brmsfit (package brms), even more options are possible for the component argument, which are not all documented in detail here. It can be any pre-defined or arbitrary distributional parameter, like mu, ndt, kappa, etc.

References

Lai, M. H. C., Kwok, O.-m., Hsiao, Y.-Y., & Cao, Q. (2018). Finite population correction for two-level hierarchical linear models. Psychological Methods, 23(1), 94–112. doi:10.1037/met0000137

Examples

data(mtcars)
m <- lm(mpg ~ wt + cyl + vs, data = mtcars)
get_varcov(m)
#>             (Intercept)         wt        cyl         vs
#> (Intercept)   11.552774  0.1680680 -1.5843752 -4.7487893
#> wt             0.168068  0.6207461 -0.3301421 -0.2797924
#> cyl           -1.584375 -0.3301421  0.3764045  0.7257641
#> vs            -4.748789 -0.2797924  0.7257641  2.6475113

# vcov of zero-inflation component from hurdle-model
data("bioChemists", package = "pscl")
mod <- hurdle(art ~ phd + fem | ment, data = bioChemists, dist = "negbin")
get_varcov(mod, component = "zero_inflated")
#>               (Intercept)          ment
#> (Intercept)  0.0115917010 -0.0009744732
#> ment        -0.0009744732  0.0001561139

# robust vcov of, count component from hurdle-model
data("bioChemists", package = "pscl")
mod <- hurdle(art ~ phd + fem | ment, data = bioChemists, dist = "negbin")
get_varcov(
  mod,
  component = "conditional",
  vcov = "BS",
  vcov_args = list(R = 50)
)
#>                   count_(Intercept)     count_phd count_femWomen
#> count_(Intercept)       0.035208897 -0.0084544838  -0.0029448560
#> count_phd              -0.008454484  0.0026430078   0.0005519599
#> count_femWomen         -0.002944856  0.0005519599   0.0114174557