How do you get confidence intervals for interactions of variables?

Question

Suppose I am building an OLS model with the following specification:

$$y = \alpha + \beta_0x_0 + \beta_1x_1 + \beta_2x_0x_1 + \epsilon$$

The variable $x_1$ is continuous and $x_0$ is binary. When $x_0$ is true the effect on $y$ of $x_1$ is $(\beta_1 + \beta_2)x_1$, but what is the confidence interval of $\beta_1 + \beta_2$?

Answer 1

The confidence interval for $\hat\beta_1$ is: $$\hat{\beta}_1 \pm t_{n-4,1-\alpha/2}\sqrt{\hat{\text{var}}(\hat\beta_1)}$$ The confidence interval for $\hat\beta_1+\hat\beta_2$, when $x_1$ is binary (0,1), is: $$(\hat\beta_1+\hat\beta_2)\pm t_{n-4,1-\alpha/2} \sqrt{\hat{\text{var}}(\hat\beta_1)+\hat{\text{var}}(\hat\beta_2)+2\hat{\text{cov}}(\hat\beta_1,\hat\beta_2)}$$ (You could look at A. Figueiras, J. M. Domenech-Massons, and Carmen Cadarso, 'Regression models: calculating the confidence intervals of effects in the presence of interactions', Statistics in Medicine, 17, 2099-2105 (1998).)

An example in R

a) Simple confidence intervals

Download http://www.stat.columbia.edu/~gelman/arm/examples/ARM_Data.zip and extract ARM_Data/earnings/heights.dta.
Prepare the dataset:

> library(foreign)                     # to import Stata data
> earnings <- read.dta("heights.dta")
> earndf <- earnings[!is.na(earnings$earn) & earnings$earn > 0, ]
> earndf$log_earn <- log(earndf$earn)
> earndf$male <- ifelse(earndf$sex == 1, 1, 0)

The model is: $$\log(\text{earning})=\alpha + \beta_0\text{height} + \beta_1\text{male} + \beta_2\text{height}\times\text{male} + \epsilon$$ Estimate the four coefficients, extract the model matrix, and calculate degrees of freedom and coefficient covariance matrix ($\sigma^2(X^TX)^{-1}$):

> mod <- lm(log_earn ~ height + male + height:male, data=earndf)
> mod_summ <- summary(mod)
> coefs <- mod_summ$coefficients[,1]; coefs
 (Intercept)       height         male  height:male 
 8.388488373  0.017007950 -0.078586216  0.007446534 
> X <- model.matrix(mod)
> dof <- nrow(X) - ncol(X)
> coefs_var <- vcov(mod)

Now you can calculate the confidence intervals:

> halfCI <- qt(0.975, dof) * sqrt(diag(coefs_var))
> matrix(c(coefs - halfCI, coefs + halfCI), nrow=4)
             [,1]        [,2]
[1,]  6.733523317 10.04345343
[2,] -0.008588732  0.04260463
[3,] -2.546456373  2.38928394
[4,] -0.029114674  0.04400774

Indeed:

> confint(mod)
                   2.5 %      97.5 %
(Intercept)  6.733523317 10.04345343
height      -0.008588732  0.04260463
male        -2.546456373  2.38928394
height:male -0.029114674  0.04400774

b) Multiple confidence intervals

To calculate the confidence interval for coefs[2] (height) plus coef[4] (height:male):

> halfCI <- qt(0.975, dof) * sqrt(coefs_var[2,2]+coefs_var[4,4]+2*coefs_var[2,4])
> as.vector(c(coefs[2]+coefs[4]-halfCI, coefs[2]+coefs[4]+halfCI))
[1] -0.00165168  0.05056065

Andrew Gelman and Jennifer Hill (Data Analysis Using Regression and Multilevel/Hierarhical Models, §7.2, where the heights example comes from) recommend another method. They summarize inferences by simulation, which gives you greater flexibility.

> library(arm)                         # the package that accompanies the book
> simul <- sim(mod, 1000)
> height_for_men <- simul@coef[,2] + simul@coef[,4]
> quantile(height_for_men, c(0.025, 0.975))
         2.5%         97.5% 
-8.938569e-05  5.006192e-02 

i.e. $(-0.00009, 0.05)$, which is not that different from $(-0.0016, 0.05)$. Simulation results vary slightly as they depend on the random number generator 'seed'. For example:

> set.seed(123)
> simul <- sim(mod, 1000)
> height_for_men <- simul@coef[,2] + simul@coef[,4]
> quantile(height_for_men, c(0.025, 0.975))
        2.5%        97.5% 
-0.001942088  0.050513401 

Answer 2

One simple trick that avoids any computation if $x_0$ is binary is to get an equivalent model. Let $z_0= 1-x_0$. It corresponds to inverting/recoding of $x_0$. Now the equation $$ y = \mu + \gamma_0 z_0 + \gamma_1 x_1 + \gamma_2 z_0 x_1 + \delta $$ has exactly the same fit as and is in fact equivalent to your equation. But the trick is that when $x_0=1$ i.e. when $z_0=0$, the effect on $y$ of $x_1$ is $\gamma_1 x_1$, which means that $\beta_1+\beta_2 = \gamma_1$ (and we can relate all beta's to all gamma's). So the inference (p-value) and the confidence interval on $\beta_1+\beta_2$ are exactly the inference and the confidence interval on $\gamma_1$. Take your favorite statistical software and you get directly your answer !