Linear combinations of coefficients in R

2018-04-12 2 minute read

Updated: 2020-12-05

One of the few features I miss from Stata is the very-intuitive lincom command. Most of the time I recreate that functionality with survey::svycontrast. But, I always forget the syntax. This code demonstrates a minimal working example.

library(survey)

N <- 100
df <- data.frame(x1 = rnorm(N), x2 = rnorm(N))
df$y <- 1 + 3 * df$x1 - 2 * df$x2 + rnorm(N, 0, 0.5)

fit <- lm(y ~ x1 + x2, data = df)

Having fit the model, we can pass unnamed vector with the right number of coefficients to get our desired linear combination:

svycontrast(fit, c(0, 1, 1))
##          contrast     SE
## contrast   1.1028 0.0698

Or a named vector with any number of coefficients (as long as the names match). Or a list of named vectors. One “gotcha” to keep in mind: I have found that the latter syntax fails for some versions of survey and may be OS-dependent.

svycontrast(fit, list(c("x2" = 1, "x1" = 1), c("x2" = 4, "x1" = 2)))
##      contrast     SE
## [1,]   1.1028 0.0698
## [2,]  -1.6700 0.2171

If you want an easy way to pass svycontrast a data.frame of where each row is a different linear combination (and each column a different variable), below is how you do that. This is usually if you want to produce similar behavior to, say, running predict on a fitted model but A) the predict call for that model doesn’t return standard errors (as with felm or fixest) and/or B) you only want to linearly combine some of the variables (equivalent to setting all other coefficients = 0). In either case, this will work:

svycontrast_df <- function(fit, newdata) {
  # Call surveyconstrast with a data frame 
  df <- newdata
  
  # Transform data.frame into a list of its row vectors
  df_list <- as.list(as.data.frame(t(newdata))) 
  df_list <- lapply(df_list, setNames, colnames(df)) # Set all character vector names inside list
  
  # Return a named list
  setNames(as.list(as.data.frame(svycontrast(fit, df_list))), c("est", "se"))
}

newdata <- data.frame(x1 = seq(-3, 3, 1))
newdata$x2 <- newdata$x^2
svycontrast_df(fit, newdata)
## $est
## [1] -26.561671 -13.832224  -4.978334   0.000000   1.102778  -1.670001  -8.318336
## 
## $se
## [1] 0.42107169 0.19615572 0.06153634 0.00000000 0.06976139 0.21706338 0.45437814