# Measuring Standard Error of two or more coefficients combined

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

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

library(alr3)
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


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(1)
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)   7.0133778 -0.67193955 -0.534382189  0.700431114  1.72549973 -1.801643527
x1dog        -0.6719396  1.13299205  0.565017653 -0.022183489 -0.04014533  0.031569426
x1mouse      -0.5343822  0.56501765  1.131288241  0.006502645  0.01176780 -0.009253944
x2happy       0.7004311 -0.02218349  0.006502645  1.395137989  1.04334118 -0.375713718
x2sad         1.7254997 -0.04014533  0.011767798  1.043341176  1.99575847 -0.679926841
x3           -1.8016435  0.03156943 -0.009253944 -0.375713718 -0.67992684  0.534679921


So the desired estimate for $$\text{SE}_{happy,dog}$$ would be sqrt(1.133 + 1.395 + 2*-0.022), or 1.576.