3
$\begingroup$

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!

$\endgroup$
7
$\begingroup$

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.

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 
| cite | improve this answer | |
$\endgroup$
5
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.