2
$\begingroup$

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

$\endgroup$
5
$\begingroup$

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 
$\endgroup$
4
  • $\begingroup$ +1 Is there a reason you are not using varcov in your R examples? $\endgroup$ – whuber Sep 2 '20 at 20:16
  • $\begingroup$ @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 :) $\endgroup$ – Sergio Sep 2 '20 at 20:49
  • $\begingroup$ 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. $\endgroup$ – whuber Sep 3 '20 at 14:13
  • $\begingroup$ @whuber Right! I did not simplify enough :) I've edited my answer. Thanks! $\endgroup$ – Sergio Sep 3 '20 at 14:30
2
$\begingroup$

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 !

$\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.