Heteroscedasticity robust variance-covariance matrix for weighted multivariate regression

Question

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))

Answer 1

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.