Confidence interval for the difference of two fitted values from a linear regression model

Question

Assume that we have a linear regression model of the form $y=\beta_0 + f_1(x_1) + f_2(x_2) + \ldots + f_n(x_n) + \epsilon$. I have written $f(x)$ to indicate that we could model the relationship between the predictors and the dependent variables flexibly, say using polynomials or splines. For simplicity's sake, let's focus on a simpler model: $$ y=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3x_2^2 + \epsilon. $$

After fitting the model to some data, we can calculate the fitted values using the estimated coefficients: $\hat{y} = \hat{\beta_0} + \hat{\beta_1} x_1 + \hat{\beta_2} x_2 + \hat{\beta_3} x_2^2$.

Now assume that we calculate the fitted values for two specific combinations of values of $x_1$ and $x_2$. Let's say we fix $x_1$ at $90$ and let $x_2 = \{2, 5\}$. That gives us two fitted values $$ \hat{y_1}=\hat{\beta_0} + \hat{\beta_1} 90 + \hat{\beta_2} 2 + \hat{\beta_3} 2^2 $$ and $$ \hat{y_2}=\hat{\beta_0} + \hat{\beta_1} 90 + \hat{\beta_2} 5 + \hat{\beta_3} 5^2 $$

Question: What's the standard error and confidence interval for the difference of these fitted values $\hat{y_2} - \hat{y_1}$?

---

Here is a simple example in R where $\beta_0 = 1.15, \beta_1 = 0.05, \beta_2 = -0.5, \beta_3 = 0.05$ and $\epsilon\sim \mathrm{N}(0, 0.25)$:

# Reproducibility
set.seed(142857)
Simulate some data
n <- 100
x1 <- rnorm(n, 100, 15)
x2 <- runif(n, 0, 10)
y <- 1.15 + 0.05x1 - 0.5x2 + 0.05*x2^2 + rnorm(100, 0, 0.5)
dat <- data.frame(y = y, x1 = x1, x2 = x2)
Fit linear regression
mod <- lm(y~x1 + poly(x2, 2, raw = TRUE), data = dat)
summary(mod)
Fitted values
predict(mod, newdata = expand.grid(x1 = 90, x2 = c(2, 5)))
       1        2 
4.885686 4.409219

Answer 1

Taking the difference of the two predicted values gives: $$ (\hat{\beta_0} + \hat{\beta_1} 90 + \hat{\beta_2} 5 + \hat{\beta_3} 5^2) - (\hat{\beta_0} + \hat{\beta_1} 90 + \hat{\beta_2} 2 + \hat{\beta_3} 2^2) = (5 - 2)\beta_2 + (5^2 - 2^2)\beta_3 = 3\beta_2 + 21\beta_3. $$ This is a linear combination of the coefficients, for which we can use the variance-covariance matrix of the model to calculate the standard error (see this Wikipedia article and this post). Specifically, let $c$ be a column vector of scalars of the same size as the coefficients in the model. Then, $c^\intercal\beta$ is a linear combination of the coefficients. The variance of $c^\intercal\beta$ is then given by: $$ \mathrm{Var}(c^\intercal\beta) = c^\intercal\Sigma c $$ where $\Sigma$ is the variance-covariance matrix of the coefficients. Taking the square root of the variance gives the standard error.

For the specific example shown in the question, we have ($c^\intercal = (0, 0, 3, 21)$) and thus:

# Reproducibility
set.seed(142857)
Simulate some data
n <- 100
x1 <- rnorm(n, 100, 15)
x2 <- runif(n, 0, 10)
y <- 1.15 + 0.05x1 + 0.05x2^2 - 0.5*x2 + rnorm(100, 0, 0.5)
dat <- data.frame(y = y, x1 = x1, x2 = x2)
Fit linear regression
mod <- lm(y~x1 + poly(x2, 2, raw = TRUE), data = dat)
summary(mod)
Linear combination of the coefficients
a <- matrix(c(0, 0, 5 - 2, 5^2 - 2^2), ncol = 1)
Standard error of the linear combination
sqrt(t(a)%%vcov(mod)%%a)
          [,1]
[1,] 0.1003602

We can check this using the emmeans package:

library(emmeans)
contrast(emmeans(mod, "x2", at = list(x1 = 90, x2 = c(2, 5))), "revpairwise", infer = c(TRUE, TRUE))
 contrast   estimate        SE df   lower.CL   upper.CL t.ratio p.value
 5 - 2    -0.4764677 0.1003602 96 -0.6756811 -0.2772542 -4.748  <.0001 

The standard error is identical.

Answer 2

An alternative approach (I agree it is devious, bit it is also interesting) is to transform your function

$$y=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3x_2^2 + \epsilon$$

into

$$y=\beta_0 + \beta_1 x_1 + \beta_2 \frac{x_2}{3} + \beta_3(x_2-2)(x_2-5) + \epsilon$$

This is the same quadratic polynomial but now you have $\hat{y}_{x_2=5} - \hat{y}_{x_2=2} = \beta_2$ and you can directly use the standard error for the coefficient $\beta_2$.