2
$\begingroup$

Say I have the following regression model:

lm(outcome ~ a * b)

                                      Estimate Std. Error    t value                         Pr(>|t|)  CI Lower
(Intercept)                         19.9300699   1.858949 10.7211470 0.000000000000000000000004040771 16.277108
a                                   -0.9772397   2.696552 -0.3624034 0.717214903624996402697888697730 -6.276148
b                                    3.2203725   2.719877  1.1840140 0.237011114082610357955971380761 -2.124370
a:b                                  3.9193899   4.133267  0.9482547 0.343491301214760680338144993584 -4.202760

Assume a and b are dummy variables.

I can compute the cell mean and standard error for the reference class by looking at the intercept. But let's say I wanted to compute the cell mean and standard errors for individuals with "a" (and not "b").

I can compute the cell mean by simple addition: 19.9300699 - 0.9772397, but what arithmetic should I perform on the standard errors to get the correct result? Similarly, what would the standard errors be associated with individuals who have both a and b?

$\endgroup$
3
  • 1
    $\begingroup$ Although this doesn't answer the question, note that executing the command lm(outcome ~ a * b - 1) in your case will automatically display what you're looking for. $\endgroup$
    – whuber
    Sep 10 '20 at 15:25
  • $\begingroup$ Right. I have been relying on re-parametrizing the model so that the intercept spits out whatever I want, but I'd really like to know how to manually compute it from the regression output (if it's even possible). Removing the intercept, for example, won't tell me what the standard error is of people with both a and b. $\endgroup$ Sep 10 '20 at 15:42
  • $\begingroup$ I don't get that last remark. Removing the intercept effectively estimates the mean separately in all four possible combinations of the binary categories a and b. $\endgroup$
    – whuber
    Sep 10 '20 at 16:34
1
$\begingroup$

As you note, the predicted cell means can be calculated by adding and subtracting the regression estimates. More precisely, in your example you're calculating a new regression parameter, $\mu_{a,\lnot b} = \beta_0 + \beta_{a}$, where $\beta_0$ is the intercept and $\beta_a$ is the slope for $a$.

The variance of this new parameter depends on the variance of $\beta_0$, the variance of $\beta_a$, and their covariance. The rule for addition is $\text{var}(x + y) = \text{var}(x) + \text{var}(y) + 2\text{cov}(x, y)$. (see Wikipedia, or my blog). The standard error of the parameter is just the square root of the variance.

Putting this together, the standard error for group a is

$$\sqrt{\text{var}(\beta_0 + \beta_a)}\\ = \sqrt{\text{var}(\beta_0) + \text{var}(\beta_a) + 2\text{cov}(\beta_0, \beta_a)}$$

Practical steps

You can find these variance and covariances from the model's covariance matrix, $\Sigma$ (although I'm not going to explain in detail what that is, except that entry Sigma[i,i] is the variance for the ith variance and Sigma[i,j] is the covariance between the ith and jth).

Finally, in the code below I go through the steps of doing these calculations manually for a model similar to your example, and show how this can be done automatically using the multcomp package.


library(tidyverse)
# Some example data
data = ToothGrowth %>%
  filter(dose >= 1) %>%
  mutate(dose = factor(dose))
data %>%
  group_by(supp, dose) %>%
  summarise(mean = mean(len), 
            sd   = sd(len),
            sem  = sd/sqrt(n()))
# A tibble: 4 x 5
# Groups:   supp [2]
# supp  dose   mean    sd   sem
# <fct> <fct> <dbl> <dbl> <dbl>
# 1 OJ    1      22.7  3.91 1.24 
# 2 OJ    2      26.1  2.66 0.840
# 3 VC    1      16.8  2.52 0.795
# 4 VC    2      26.1  4.80 1.52 

contrasts(data$supp) = c(0, 1)
contrasts(data$dose) = c(0, 1)
m = lm(len ~ supp * dose, data=data)
summary(m) %>% coef() %>% round(2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)    22.70       1.14   19.97     0.00
# supp1          -5.93       1.61   -3.69     0.00
# dose1           3.36       1.61    2.09     0.04
# supp1:dose1     6.01       2.27    2.64     0.01

b = coef(m)
sigma = vcov(m)
#             (Intercept) supp1 dose1 supp1:dose1
# (Intercept)        1.29 -1.29 -1.29        1.29
# supp1             -1.29  2.58  1.29       -2.58
# dose1             -1.29  1.29  2.58       -2.58
# supp1:dose1        1.29 -2.58 -2.58        5.17

# Let's get supp==OJ, dose==2
# = intercept + b_dose
estimate = b[1] + b[2] # 16.77
sem = sqrt(sigma[1,1] + sigma[2,2] + 2*sigma[1,2]) # 1.136771

# You can do this automatically using the multcomp package
library(multcomp)
ht = glht(m, rbind(c(1, 1, 0, 0)))
summary(ht)

#    Simultaneous Tests for General Linear Hypotheses
# Fit: lm(formula = len ~ supp * dose, data = data)
# 
# Linear Hypotheses:
#        Estimate Std. Error t value Pr(>|t|)    
# 1 == 0   16.770      1.137   14.75   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)

# Alternative code
indices = c(1, 2)
estimate = sum(b[indices]) # 16.77
v = 0
for(i in indices){
  for(j in indices){
    if(i == j) v = v + sigma[i, j]
    if(i > j) v = v + 2*sigma[i, j]
  }
}
sem = sqrt(v) # 1.13

# Now let's do supp==CV, dose==2
# = intercept + b_supp + b_dose + interaction
indices = c(1,2,3,4)
estimate = sum(b[indices]) # 26.14
v = 0
for(i in indices){
  for(j in indices){
    if(i == j) v = v + sigma[i, j]
    if(i > j) v = v + 2*sigma[i, j]
  }
}
sem = sqrt(v) # 1.13

# Using multcomp
glht(m, rbind(c(1, 1, 1, 1))) %>% summary()
#    Simultaneous Tests for General Linear Hypotheses
# Fit: lm(formula = len ~ supp * dose, data = data)
# 
# Linear Hypotheses:
#        Estimate Std. Error t value Pr(>|t|)    
# 1 == 0   26.140      1.137      23   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)
```
$\endgroup$
3
  • 1
    $\begingroup$ After writing this, it occurs to me that this question has almost certainly already been asked and answered on this site already. $\endgroup$
    – Eoin
    Sep 10 '20 at 16:30
  • 1
    $\begingroup$ Here's one search that will pick up a few of those answers: stats.stackexchange.com/…. $\endgroup$
    – whuber
    Sep 10 '20 at 16:35
  • $\begingroup$ This is perfect. Thank you so much. $\endgroup$ Sep 10 '20 at 16:38

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.