Computes the asymptotic HC0 covariance matrix for multiple models fit to the same data, with cross-model covariances included. These are equivalent to the covariance computed when stacking models using M-estimation.

vcovSUEST.R

# Computes joint HC0 covariance matrix of several models fit to the same data.
# `fits` should be a list of model fits (e.g., output of a call to lm or glm, etc.)
# To include models fit to subsets of data, fit models to whole dataset with weights
# close to 0 for units to be excluded. Relies on `sandwich` functionality. Returns
# a symmetric matrix with no dimnames. Individual model covariances are on the block
# diagonals; between-model covariances are on the off-diagonals. See
# https://github.com/kylebutts/vcovSUR for a more mature implementation. See Mize
# et al. (2019) <https://doi.org/10.1177/0081175019852763> for theory and
# application.

vcovSUEST <- function(fits) {
    
    inf_func <- lapply(fits, function(f) { 
        b <- sandwich::bread(f)
        ef <- sandwich::estfun(f)
        inf <- tcrossprod(b, ef)/nobs(f)
        inf[is.na(inf)] <- 0
        inf
    })
    
    coef_lengths <- unlist(lapply(inf_func, ncol))
    
    l <- c(0, cumsum(coef_lengths))
    
    coef_inds <- lapply(seq_along(fits), function(i) seq_len(coef_lengths[i]) + l[i])
    
    #VCOV matrix to be returned
    V <- matrix(NA_real_, nrow = sum(coef_lengths), ncol = sum(coef_lengths))
    
    for (i in seq_along(fits)) {
        ind_i <- coef_inds[[i]]
        
        #Usual within-model HC0 vcov
        V[ind_i, ind_i] <- tcrossprod(inf_func[[i]])
        
        for (j in seq_along(fits)[-seq_len(i)]) {
            ind_j <- coef_inds[[j]]
            
            #between-model vcov components
            V[ind_i, ind_j] <- tcrossprod(inf_func[[i]], inf_func[[j]])
            V[ind_j, ind_i] <- t(V[ind_i, ind_j])
        }
    }
    
    V
}