Covariance matrix of the fixed-effect estimates

Covers brms's population-level coefficients and names its rows as brms does (Intercept, sigma_Intercept, x). A distributional parameter nobody wrote a formula for is not one of them: frmtmb estimates sigma or nu as an intercept-only linear predictor, brms reports it as a parameter of its own, and summary() shows it under Further Distributional Parameters. Dpars fixed to constants are excluded too.

Usage

# S3 method for class 'frmtmb_fit'
vcov(
  object,
  correlation = FALSE,
  pars = NULL,
  ...,
  full = FALSE,
  cluster = NULL,
  type = "CR0"
)

Arguments

object: A frmtmb_fit.
correlation: If TRUE, the correlation matrix instead, as in brms.
pars: Row and column names to keep, in the order given, as in brms.
...: Refused: an argument the method does not have is an error naming it, rather than silently changing nothing.
full: If TRUE, include covariance parameters (theta), named as in confint() (the glmmTMB vcov(full = TRUE) convention). full = TRUE keeps the INTERNAL parameters, because it is the matrix a delta-method calculation on confint()'s rows needs; the default block takes brms's names. An ordinal fit's thresholds are named there as confint() names them, by what they are: Intercept[1], delta, or a transform such as log(Intercept[2] - Intercept[1]).
cluster: Optional clustering factor. When given, the result is vcov_cluster()'s cluster-robust covariance instead of the model-based one.
type: Small-sample correction for cluster, see vcov_cluster().

Value

A covariance matrix.

Details

full = TRUE is the joint covariance of the whole outer parameter vector on its internal scale: the fixed-effect coefficients, the covariance parameters theta (log standard deviations, Fisher-z correlations, and whatever else a structure keeps there), and any extra parameters such as the ordinal thresholds. It is the matrix a delta-method calculation on a variance component needs, and it is what hypothesis() uses for method = "wald" - so an ICC or a heritability is usually easier to ask for through hypothesis(), which names the components for you, than to assemble by hand from this matrix.

Under REML = TRUE (or frmtmb_control(profile = TRUE)) the fixed effects are integrated out of the outer problem, so they are not part of full = TRUE; the block comes from the joint precision and carries exactly the parameters confint.frmtmb_fit() reports. vcov(object) is still the fixed-effect covariance there.

Under REML that covariance is the beta block of the inverse joint precision. It is the conditional covariance at the estimated variance parameters, sigma^2 (X' V^-1 X)^-1, plus J V_theta J': the variance the uncertainty of the variance parameters theta adds, with J the derivative of the estimate of beta in theta and V_theta their covariance. That is glmmTMB's definition, and its REML standard errors agree with frmtmb's to ten digits. lme4's vcov() reports the conditional part alone, so the two differ where J is not zero: on an unbalanced design by up to 1e-3 of the standard error, on a balanced one not at all. On an ill-conditioned design both parts are taken in the design's coordinates and held as a factor; the products warning of an ML fit applies.

Passing cluster forwards to vcov_cluster() for the cluster-robust (sandwich) covariance, in the sandwich::vcovCL() spelling: vcov(fit, cluster = ~ g, type = "CR1").

Examples

set.seed(1)
dd <- data.frame(x = rnorm(100), g = factor(rep(1:10, 10)))
dd$y <- rnorm(100, 1 + 0.5 * dd$x + rnorm(10, 0, 0.8)[dd$g], 1)
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd)

# standard errors of the fixed effects, which is fixef()'s Est.Error
sqrt(diag(vcov(fit)))
#> Intercept         x 
#> 0.2745684 0.1132853 
# the covariance parameters join the block on their internal scale,
# under the internal names, which is confint()'s vocabulary
rownames(vcov(fit, full = TRUE))
#> [1] "(Intercept)"       "x"                 "sigma_(Intercept)"
#> [4] "theta_1"          

# the matrix is what a delta-method calculation needs
V <- vcov(fit)
a <- c(1, 2)                       # prediction at x = 2, no group
sqrt(drop(t(a) %*% V[1:2, 1:2] %*% a))
#> [1] 0.3480407