Confidence interval over point estimate given regression parameters

Question

In Bayesian analysis, the posterior distribution is often sampled when the PDF is not intractable (often). If the samples are of length $n$, then every index in $range(1,n)$ corresponds to a valid sample from that (often multi-dimensional) distribution of simple linear regression; the slope, intercept and noise are the three dimensions of this posterior distribution (assuming no hyper-priors.)

To compute the credible distribution around a given point $P(\hat{y}|x, \beta_0, \beta_1)$, a simulation can be run where indices from $range(1,n)$ can be sampled and the parameter values retrieved; subsequently, the value of $\hat{y}$ can be computed and appended to an array. Lastly, the array could be sorted and processed as desired (ex 95% credible interval.)

My question is, what's the analogous process for Frequentist statistics? Even if the slope or coefficient is statistically significant and the confidence intervals for each are known, this isn't the joint distribution across model parameters. So, I suspect that one cannot sample parameters from each of the confidence intervals independently to compute the distribution over the point estimate.

Is this possible from the Frequentist toolkit? Or is asking such a question a purely Bayesian construct, and Frequentists would not be interested beyond knowing what the parameter values mean and whether they are significant?

Answer 1

The variance-covariance matrix for the coefficients contains the covariances in the off-diagonal cells. This allows you to construct a linear combination of random model coefficients and the covariates and its variance, which in turn gets you the confidence/coverage interval for the expected value or a prediction interval.

If you are interested in other parts of the distribution, you can use quantile regression to do something similar.

---

Example in R

# Load necessary library
library(stats)

# Create some synthetic data
set.seed(123) # for reproducibility
x <- 1:100
y <- 2 + 3 * x + rnorm(100, mean = 0, sd = 20) # y = 2 + 3x + noise

# Fit a linear regression model
model <- lm(y ~ x)

# Display the summary of the model
summary(model)

# Extract the variance-covariance matrix of the model coefficients
vcov_matrix <- vcov(model)

# Print the variance-covariance matrix
print(vcov_matrix)

# Calculated expected y for x = 2 with 95% confidence interval
new_data <- data.frame(x = 2)
prediction <- predict(model, new_data, interval = "confidence", level = 0.95)

# Print the expected value with with 95% confidence interval
print(prediction)

You could also get the predictional interval for the outcome, rather than for its expected value with

prediction <- predict(model, new_data, interval = "prediction", level = 0.95)

The prediction will be the same as above, but the interval will be wider since the outcome has a higher variance than the mean outcome. More on this distinction here.