Measuring Standard Error of two or more coefficients combined

Question

A silly example.

x1 <- as.factor(c(rep("dog", 3), rep("cat", 3), rep("mouse", 3)))
x2 <- as.factor(rep(c("happy", "sad", "angry"), 3))
x3 <- rnorm(9, 0, 1) + runif(9, 3, 5)
y <- rnorm(9, 10, 2)

Call:
lm(formula = y ~ x1 + x2 + x3)

Residuals:
       1        2        3        4        5        6        7        8        9 
-1.57949  1.51090  0.06859  1.59378 -2.50472  0.91094 -0.01429  0.99383 -0.97953 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 19.69066    6.98359   2.820   0.0668 .
x1dog       -5.89792    3.40777  -1.731   0.1819  
x1mouse     -2.05016    3.32847  -0.616   0.5815  
x2happy      1.57757    1.99707   0.790   0.4872  
x2sad       -0.02729    2.09737  -0.013   0.9904  
x3          -1.83281    1.19873  -1.529   0.2237  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.336 on 3 degrees of freedom
Multiple R-squared:  0.6874,    Adjusted R-squared:  0.1664 
F-statistic: 1.319 on 5 and 3 DF,  p-value: 0.4362

Say I wanted to know the marginal effect of a 'happy dog' on y. I would add x1dog and x2happy, and say something like, "compared to an angry cat, a happy dog has a -4.32 marginal effect on y."

My question is, what standard error would I give on this estimate? I don't think I should just add the two respective SE's. What is the approach?

Thanks!

Answer 1

1) A linear combination a of the coefficients has a standard error of $\sqrt(a'Va)$ where $V$ is the variance covariance matrix of the coefficients. (See this wikipedia link on the variance of a linear combination of random variables.) We can get $V$ using vcov in R so using fm from the reproducible code in the Note at the end we have:

a <- c(0, 1, 0, 1, 0, 0)
c(sqrt(t(a) %*% vcov(fm) %*% a))
## [1] 2.940084

2) delta method The delta method can be used to get an approximate standard error of a general (differentiable) nonlinear function of the coefficients. When the function is linear, as it is here, then the Taylor series approximation underlying the delta method is exact and gives the same result as above.

# update - alr3 no longer on CRAN but same code works with car
# library(alr3)
library(car)
deltaMethod(fm, "x1dog+x2happy")
##                   Estimate       SE     2.5 %   97.5 %
## x1dog + x2happy -0.2095575 2.940084 -5.972016 5.552901

3) general linear hypothesis Testing the hypothesis $a'\beta = 0$ (where $\beta$ is the coefficient vector) involves using the desired standard error in the denominator of the t-test suggesting that we look for a function to perform such hypothesis test. In fact, running glht with the hypothesis x1dog + x2happy = 0 does display the desired standard error, as shown below:

library(multcomp)
summary(glht(fm, "x1dog + x2happy = 0"))

giving the following where the Std. Error shown is the standard error of x1dog + x2 happy. It equals the standard error computed in (1) and (2).

         Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = y ~ x1 + x2 + x3)

Linear Hypotheses:
                     Estimate Std. Error t value Pr(>|t|)
x1dog + x2happy == 0  -0.2096     2.9401  -0.071    0.948
(Adjusted p values reported -- single-step method)

Note

We assume the input is the following where we use set.seed to make it reproducible.

set.seed(123)
x1 <- as.factor(c(rep("dog", 3), rep("cat", 3), rep("mouse", 3)))
x2 <- as.factor(rep(c("happy", "sad", "angry"), 3))
x3 <- rnorm(9, 0, 1) + runif(9, 3, 5)
y <- rnorm(9, 10, 2)
fm <- lm(y ~ x1 + x2 + x3)

coef(fm)
## (Intercept)       x1dog     x1mouse     x2happy       x2sad          x3 
##  12.4103736   0.1637373  -1.3375086  -0.3732948  -1.4833220  -0.3470490 

Answer 2

The sum of two correlated (jointly-)normal variables, $Z=X+Y$, is itself normal, $N(\mu_{X+Y},\sigma_{X+Y})$.

The expectation of the distribution is the simple sum, $\mu_{X+Y} = \mu_X+\mu_Y$.

The S.D. of the distribution is not quite as simple since you need to know the covariance of the two distributions, but still straightforward to calculate: $\sigma_{X+Y}=\sqrt{\sigma_X^2 + \sigma_Y^2 + 2*\sigma_{XY}}$.

Now all we need is the covariance $\sigma_{XY}$, which is available in R using vcov() (for variance-covariance matrix) on most model objects.

set.seed(123)
x1 <- as.factor(c(rep("dog", 3), rep("cat", 3), rep("mouse", 3)))
x2 <- as.factor(rep(c("happy", "sad", "angry"), 3))
x3 <- rnorm(9, 0, 1) + runif(9, 3, 5)
y <- rnorm(9, 10, 2)


fit <- lm(formula = y ~ x1 + x2 + x3)

covariances <- vcov(fit)

            (Intercept)      x1dog   x1mouse    x2happy      x2sad         x3
(Intercept)   29.605163 -3.6857520 -7.793895 -7.6087847 -7.7306824 -4.4231039
x1dog         -3.685752  3.2430017  2.033472  0.4826921  0.4923812  0.3515720
x1mouse       -7.793895  2.0334723  4.520103  1.4051380  1.4333433  1.0234414
x2happy       -7.608785  0.4826921  1.405138  4.4357063  2.9270107  0.9931674
x2sad         -7.730682  0.4923812  1.433343  2.9270107  4.4909978  1.0131032
x3            -4.423104  0.3515720  1.023441  0.9931674  1.0131032  0.7233800

So the desired estimate for $\text{SE}_{happy,dog}$ would be sqrt(3.243... + 4.435... + 2*0.482...), or 2.940084.

Update

Changed sample data to show that this method matches results from the functions recommended by the accepted answer.