Function to compute simple linear regression coefficients and their covariances

coef_cov.R

#####
# Function: coef_cov
# Author: Brenton M. Wiernik
# Date: 2020-01-27
# Version: 1.0
# 
# Computes univariate linear regresssion models,
# predicting each response variable [`.y`] with each
# predictor variable [`.x`]. Returns the coefficients
# and their variance-covariance matrix.
# 
# Arguments:
# .data  A data frame containing variables to estimate models. 
# .y     The column names or indices in .data for the response variables.
# .x     The column names or indices in .data for the predictor variables.
# 
# Returns: 
# A list containing the estimated coefficients (`coef`) and their 
# variance-covariance matrix (`vcov`)
# 
# Note: Missing values in `.data` will be dropped listwise.
# 
#####

coef_cov <- function(.data, .y, .x) {
  if (is.null(colnames(.data))) stop(".data must have names")
  if (is.numeric(.y)) y <- colnames(.data[,.y]) else y <- .y
  if (is.numeric(.x)) x <- colnames(.data[,.x]) else x <- .x
  dat <- na.omit(.data[,c(y, x)])
  m <- sapply(x, function(.) lm(as.matrix(dat[,y]) ~ dat[,.]), simplify = FALSE)
  COEF <- as.vector(sapply(m, function(.) coef(.)[seq(2, 2*length(y), by = 2)]))
  names(COEF) <- as.vector(outer(y, x, "paste", sep="_"))
  vcs <- sapply(m, function(.) vcov(.)[seq(2, 2*length(y), by = 2), seq(2, 2*length(y), by = 2)], simplify = FALSE)
  eps <- sapply(m, resid, simplify = FALSE)
  VC <- as.matrix(Matrix::bdiag(vcs))
  rownames(VC) <- colnames(VC) <- as.vector(outer(y, x, "paste", sep="_"))
  VC[upper.tri(VC)] <- 0
  for (r in seq(1, nrow(VC), length(y))[-1]) {
    rs <- r:(r + length(y) - 1)
    r_x <- which(r == seq(1, nrow(VC), length(y)))
    for (c in seq(1, ncol(VC), length(y))[-length(x)]) {
      cs <- c:(c + length(y) - 1)
      c_x <- which(c == seq(1, nrow(VC), length(y)))
      xr <- dat[,x[r_x]]
      xc <- dat[,x[c_x]]
      VC[rs, cs] <- solve(t(xr) %*% xr) %*% t(xr) %*% eps[[r_x]] %*% t(eps[[c_x]]) %*% xc %*% solve(t(xc) %*% xc)
    }
  }
  VC <- VC + t(VC) - diag(diag(VC))
  return(list(coef = COEF, vcov = VC))
}

coef_cov(psych::bfi, 1:5, 6:11)