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
mod
via thevcov
function. $\endgroup$