# Heteroscedasticity robust variance-covariance matrix for weighted multivariate regression

I need heteroscedasticity robust standard errors for a multivariate linear model (MLM) with weights. In R we usually use sandwich::vcovHC with type "HC0" to get the White variance-covariance matrix. But while this works for univariate models with weights or MLMs without weights, it doesn't seem to support MLMs with weights.

## Example

fit1 <- lm(cbind(X1, X2) ~ . - w, dat)
sandwich::vcovHC(fit1, type='HC0')  ## unweighted MLM, works

fit2 <- lm(cbind(X1, X2) ~ . - w, dat, weights=w)
sandwich::vcovHC(fit2, type='HC0') ## weighted MLM, fails
# Error in SSD.mlm(object) : 'mlm' objects with weights are not supported


However, since it does work for the weighted univariate case,

fit31 <- lm(X1 ~ . - X2 - w, dat, weights=w)
lmtest::coeftest(fit31, sandwich::vcovHC(fit31, type='HC0'))
fit32 <- lm(X2 ~ . - X1 - w, dat, weights=w)
lmtest::coeftest(fit32, sandwich::vcovHC(fit32, type='HC0'))


there should be an easy way to implement this for MLMs.

In the source code, the error boils down to computing the "matrix of residual sums of squares and products."

stats:::SSD.mlm(fit2)
# Error in stats:::SSD.mlm(fit2) :
#   'mlm' objects with weights are not supported


## Statistics

stats:::SSD.mlm essentially calculates crossprod(object$residuals)―the crossproduct of the residuals, $$r = Y - \hat{Y}\quad (1)$$ $$SSD=r'r \quad(2)$$ but it is obviously programmed to stop if there are weights in the MLM object. I don't understand why sandwich:::vcovHC.mlm was not implemented for MLMs with weights, since it was possible for MLMs without weights. Do I assume correctly that it would not be sufficient to simply adapt equation (2) for the case with weights, because that would have been too easy to implement? So my question is, how can we adapt the code, or which equations do we need to adapt the code for MLMs? A solution would be very valuable because we could compute robust errors for tens of thousands of regressions at once instead of running a loop for hours. Data: set.seed(42); n <- 2 + 2; m <- 10 dat <- data.frame(matrix(rnorm(m*n), m, n)) dat$w <- runif(nrow(dat))


You can at least get some type of result if you use the data transformation that's equivalent to weighted least squares when it's assumed that there are no correlations among error terms. As this answer suggests:

If you pre-multiply each of the design matrix and the outcome vector by the diagonal matrix of the square roots of the case weights, then OLS gives the same result as weighted least squares.

Starting with your data frame dat and the unweighted model fit1, get the diagonal matrix containing the square roots of the weights, then pre-multiply each of the outcome matrix and the unweighted model matrix by that diagonal matrix.

dsw <- diag(sqrt(dat$w)) wtOut <- dsw %*% cbind(dat$X1,dat\$X2)
colnames(wtOut) <- c("X1","X2")
wtMM <- dsw %*% model.matrix(fit1)


As the weighted model matrix already contains the appropriate coding for the intercept, now just do lm() without an additional intercept.

wtMLM <- lm(wtOut ~ wtMM -1)
wtMLM
#
# Call:
# lm(formula = wtOut ~ wtMM - 1)
#
# Coefficients:
#                  X1        X2
# wtMM(Intercept)   0.50188   0.28926
# wtMMX3            0.09423  -0.34577
# wtMMX4            0.01773   0.18896


Those coefficient estimates are identical to those for your weighted fit2 (not shown). As lm() didn't know that there were weights involved, you can get either a standard vcov() or a sandwich estimate.

dim(vcov(wtMLM))
# [1] 6 6
dim(sandwich::vcovHC(wtMLM, type='HC0'))
# [1] 6 6


I nevertheless worry, as in the linked answer above, whether this type of weighting makes sense in the multivariate regression context.

• Thanks, satisfying how closely coefficients and statistics agree with univariate WLS. I first did something very similar using this answer that calculates WLS per hand. However, there is no obvious square root involved, where in your method there is, could you explain that? Second, can you elaborate on the assumption that the errors are not correlated? The strength of multivariate models is that it also considers the correlation of the covariates across outcomes, does this still apply here? Jul 3, 2023 at 8:09
• @jay.sf the Wikipedia page shows the equivalence of usual weighting and the square-root factorization approaches. The manual approach you cite would have to be extended for multivariate outcomes. Multivariate models consider across-outcome correlations, but the weight matrix must be diagonal for the simple square-root factorization to work. It assumes that within-outcome correlations, as handled by generalized least squares, aren't at work. The sandwich vcov can help with that.
– EdM
Jul 3, 2023 at 12:22