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

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.05*x1 - 0.5*x2 + 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

• Just do the usual covariance calculation: after all, the difference is a linear combination of $\hat\beta$ and you can extract the variance-covariance matrix of that difference from mod via the vcov function. – whuber Oct 15 '20 at 20:44

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.05*x1 + 0.05*x2^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.

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