# How do you get confidence intervals for interactions of variables?

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$$?

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

Prepare the dataset:

> library(foreign)                     # to import Stata data
> 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

• +1 Is there a reason you are not using varcov in your R examples? – whuber Sep 2 '20 at 20:16
• @whuber Which varcov? However, in a first draft I used... nothing, e.g. RMS <- (t(y) %*% (diag(nrow(X)) - X %*% solve(t(X) %*% X) %*% t(X)) %*% y) / (nrow(X)-ncol(X)), then I've tried to simplify my code :) – Sergio Sep 2 '20 at 20:49
• Sorry, I meant vcov, which is part of base R. It will be numerically more stable than your solution but, more importantly, its use would greatly clarify your calculations. – whuber Sep 3 '20 at 14:13
• @whuber Right! I did not simplify enough :) I've edited my answer. Thanks! – Sergio Sep 3 '20 at 14:30

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 !