3
$\begingroup$

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
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$ – whuber Oct 15 '20 at 20:44
3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

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