Skip to main content
Cross Validated

Return to Answer

Added missing square
Source Link
Martin Modrák
Martin Modrák
  • 3.1k
  • 14
  • 30

We can realize that if $\sigma$ was known, the distribution of the MLE of the slope is $\frac{\sigma}{\sum (x_i - m)^2}$$\frac{\sigma^2}{\sum (x_i - m)^2}$ where $m$ is the observed mean of $x$ (taken from Wiki )

So the prediction interval will be derived from t distribution centered at the MLE from the first batch of data and with standard error of $\sqrt{s^2 + \frac{\hat\sigma}{\sum (\bar{x}_i - \bar{m})^2}}$$\sqrt{s^2 + \frac{\hat{\sigma}^2}{\sum (\bar{x}_i - \bar{m})^2}}$ where $\hat\sigma$ is the estimate of $\sigma$ from the first regression, $s$ is the standard error from the first regression and $\bar{m}$ is the mean of $\bar{x}$.

We can realize that if $\sigma$ was known, the distribution of the MLE of the slope is $\frac{\sigma}{\sum (x_i - m)^2}$ where $m$ is the observed mean of $x$ (taken from Wiki )

So the prediction interval will be derived from t distribution centered at the MLE from the first batch of data and with standard error of $\sqrt{s^2 + \frac{\hat\sigma}{\sum (\bar{x}_i - \bar{m})^2}}$ where $\hat\sigma$ is the estimate of $\sigma$ from the first regression, $s$ is the standard error from the first regression and $\bar{m}$ is the mean of $\bar{x}$.

We can realize that if $\sigma$ was known, the distribution of the MLE of the slope is $\frac{\sigma^2}{\sum (x_i - m)^2}$ where $m$ is the observed mean of $x$ (taken from Wiki )

So the prediction interval will be derived from t distribution centered at the MLE from the first batch of data and with standard error of $\sqrt{s^2 + \frac{\hat{\sigma}^2}{\sum (\bar{x}_i - \bar{m})^2}}$ where $\hat\sigma$ is the estimate of $\sigma$ from the first regression, $s$ is the standard error from the first regression and $\bar{m}$ is the mean of $\bar{x}$.

Source Link
Martin Modrák
Martin Modrák
  • 3.1k
  • 14
  • 30

Thanks to Jarle Tufto who provided an outline of the main idea in the comments. All errors are my own.

We can realize that if $\sigma$ was known, the distribution of the MLE of the slope is $\frac{\sigma}{\sum (x_i - m)^2}$ where $m$ is the observed mean of $x$ (taken from Wiki )

This means that the difference of the two MLEs is also going to be normally distributed. We however do not know $\sigma$, and plugging in the estimate of sigma from the first part, the resulting quantity is $t$-distributed with $n−2$ degrees of freedom, where $n$ is the sample size of the observed data.

So the prediction interval will be derived from t distribution centered at the MLE from the first batch of data and with standard error of $\sqrt{s^2 + \frac{\hat\sigma}{\sum (\bar{x}_i - \bar{m})^2}}$ where $\hat\sigma$ is the estimate of $\sigma$ from the first regression, $s$ is the standard error from the first regression and $\bar{m}$ is the mean of $\bar{x}$.

Adapting the code example from the question:



x <- rep(0:5, times = 4)
xbar <- rep(0:5, times = 2)

N_sims <- 10000

#Pre-allocate
y <- matrix(NA_real_, nrow = N_sims, ncol = length(x))
ybar <- matrix(NA_real_, nrow = N_sims, ncol = length(xbar))
b1 <- numeric(N_sims)
b1ci <- matrix(NA_real_, nrow = N_sims, ncol = 2)
b1bar <- numeric(N_sims)

b1pred <- matrix(NA_real_, nrow = N_sims, ncol = 2)

for(i in 1:N_sims) {
  # Choose the true values randomly
  # The results are quite similar whether I choose
  # those separately for each simulation or
  # once at the begining of the script
  true_b0 <- 1 + rnorm(1)
  true_b1 <- 0.5 + rnorm(1)
  true_sigma <- 1 + rlnorm(1)


  y_ <- rnorm(length(x), true_b0 + x * true_b1, sd = true_sigma)
  y[i,] <- y_
  fit <- lm(y_ ~ x)
  b1[i] <- coef(fit)["x"]
  b1ci[i,] <- confint(fit)["x",]

  # Calculate the prediction interval for b1bar
  b1bar_var <- sigma(fit)^2 / sum((xbar - mean(xbar))^2)
  b1_var <- vcov(fit)["x", "x"]
  diff_se <- sqrt(b1bar_var + b1_var)
  b1pred[i, ] <- b1[i] + qt(c(0.025,0.975), df = length(x) - 2) * diff_se

  ## Actually fit the new data
  ybar_ <- rnorm(length(xbar), true_b0 + xbar * true_b1, sd = true_sigma)
  ybar[i,] <- ybar_
  fitbar <- lm(ybar_ ~ xbar)
  b1bar[i] <- coef(fitbar)["xbar"]
}

cat("Coverage prediction: ", mean(b1bar >= b1pred[,1] & b1bar <= b1pred[,2]))

Shows the expected coverage behaviour.