# What is the Matrix form of the prediction interval for Weighted Least Squares Regression

For ordinary least squares, I know the confidence interval is $$\hat{Y} ± t\cdot\sqrt{\hat{\sigma}^2 * x_{new}(X^{\top}X)^{-1}x_{new}^{\top}}$$ and $$X$$ is the design matrix, n is the number of observations, and p is the number of parameters in the model, and:

$$\hat{\sigma}^2 = \frac{Y^{\top} \left(I-X(X^{\top}X)^{-1}X^{\top} \right)Y}{n-p}$$

My design Matrix has observations in rows.

The prediction interval is similarly defined as: $$\hat{Y} ± t\cdot\hat{\sigma} \cdot \sqrt{ x_{new}(X^{\top}X)^{-1}x_{new}^{\top} + 1/q}$$ where q is the number of replicate measurements of the future sample.

For weighted least squares, the confidence interval is: $$\hat{Y} ± t\cdot\sqrt{\hat{\sigma}^2\cdot x_{new}(X^{\top}WX)^{-1}x_{new}^{\top}}$$

$$\hat{\sigma}^2 = \frac{Y^{\top}W\left(I-X(X^{\top}WX)^{-1}X^{\top}W \right)Y}{n-p}$$

What is the equation for the weighted prediction interval. I would guess that instead of 1/q, the 1 would be dependent on the overall scale of the weight matrix in some way, or the weight matrix should be scaled such that the squared sum of the elements is one, or something like that.

Doing some tests, the value of the confidence interval does not appear to change if the weight matrix is multiplied by a scalar, but If the prediction interval is calculated as: $$\hat{Y} ± t*\hat{\sigma} \sqrt{ x_{new}(X^{\top}WX)^{-1}x_{new}^{\top} + 1/q}$$, then scaling the weights does affect the interval.

• Could you explain the $1/q$ term or give a source where you got this from. The typical expression is $$\hat{Y} \pm t\cdot \hat{\sigma} \sqrt{1/n +1}$$ where $n$ refers to the number of past observations instead of the number of future observations. The $1/n$ relates to the error of the estimate in $\hat{Y}$ which scales with $\sqrt{1/n}$ (or in the case of OLS with $\sqrt{x_n(X^TX)^{-1}x_n^T}$ instead of $1/n$). en.m.wikipedia.org/wiki/… Commented Aug 14, 2022 at 11:08
• @SextusEmpiricus I believe the 1/q term was included in the equation of the prediction interval under simple linear regression in the textbook Quantitative Chemical Analysis by Daniel Harris Commented Aug 15, 2022 at 13:56
• @SextusEmpiricus I've modified the question slightly. The text refers to it as the number of replicate measurements of the future sample. Commented Aug 15, 2022 at 14:21
• @KevinNowaczyk $X(X^\top WX)^{-1}X^\top W$ is not symmetric. But since $-Y^\top W^\top X(X^\top WX)^{-1}X^\top WY = (-2+1)[Y^\top W^\top X(X^\top WX)^{-1}X^\top WY]$ which in turn equals $-2\times\hat\beta^\top X^\top WY+1\times\hat\beta^\top X^\top WX\hat\beta$, your expression for the estimator of $\sigma^2$ in the weighted least squares case is in fact equivalent to the one in my answer. Commented Aug 15, 2022 at 14:50
• @KevinNowaczyk Are you referring to an equation given in the section on Calibration Curves? Commented Aug 17, 2022 at 10:21

Trying to match your notation and what you have written already, I'll assume that $$Y=X\beta+\varepsilon,\varepsilon\sim\mathcal{N}\left(0,\sigma^{2}W^{-1}\right)$$ for a known positive-definite and diagonal matrix $$W$$ and that $$\hat{\beta}$$ is the weighted least squares estimator $$\left(X^{\top}WX\right)^{-1}X^{\top}WY$$.
This gives \begin{align*} \hat{\beta} & \sim\mathcal{N}\left(\beta,\sigma^{2}\left(X^{\top}WX\right)^{-1}\right),\\ x_{\mathrm{new}}\hat{\beta} & \sim\mathcal{N}\left(x_{\mathrm{new}}\beta,\sigma^{2}x_{\mathrm{new}}\left(X^{\top}WX\right)^{-1}x_{\mathrm{new}}^{\top}\right),\\ Y_{\mathrm{new}}-x_{\mathrm{new}}\hat{\beta} & \sim\mathcal{N}\left(0,\sigma^{2}w_{\mathrm{new}}^{-1}+\sigma^{2}x_{\mathrm{new}}\left(X^{\top}WX\right)^{-1}x_{\mathrm{new}}^{\top}\right), \end{align*} where $$w_{\mathrm{new}}^{-1}$$ is the diagonal entry of $$W^{-1}$$ corrsponding to the future observation $$Y_{\mathrm{new}}$$.
Standardizing and replacing $$\sigma^{2}$$ with the unbiased estimator $$\hat{\sigma}^{2}=\left(Y-X\hat{\beta}\right)^{\top}W\left(Y-X\hat{\beta}\right)/\left(n-p\right)$$ yields $$\frac{Y_{\mathrm{new}}-x_{\mathrm{new}}\hat{\beta}}{\hat{\sigma}\sqrt{w_{\mathrm{new}}^{-1}+x_{\mathrm{new}}\left(X^{\top}WX\right)^{-1}x_{\mathrm{new}}^{\top}}}\sim t_{n-p}$$ and thus the prediction interval $$x_{\mathrm{new}}\hat{\beta}\mp t_{n-p}^{1-\alpha/2}\hat{\sigma}\sqrt{w_{\mathrm{new}}^{-1}+x_{\mathrm{new}}\left(X^{\top}WX\right)^{-1}x_{\mathrm{new}}^{\top}},$$ where $$t_{n-p}^{1-\alpha/2}$$ is the $$\left(1-\alpha/2\right)$$ quantile of the $$t$$-distribution with $$n-p$$ degrees of freedom.
This prediction interval contains the future observation $$Y_{\mathrm{new}}$$ with probability $$1-\alpha$$. In practice however, it can only be calculated in this way if $$w_{\mathrm{new}}^{-1}$$ has a known functional relationship with $$x_{\mathrm{new}}$$.
• +1 One sidenote: You start with a general problem statement where $W$ is a positive-definite matrix (and not more restricted where $W$ is a diagonal matrix). Then shouldn't we also consider potential correlation between the $\epsilon_{new}$ and the errors in $\hat{\beta}$? Commented Aug 14, 2022 at 22:16
• @SextusEmpiricus I'm considering the weighted least squares case (where $W$ is diagonal) as opposed to the generalized least squares (GLS) case (where $W$ doesn't have to be diagonal), it's now stated explicitly$−$sorry if this wasn't clear. In the general GLS case, independence of $\varepsilon_\mathrm{new}$ and $\varepsilon$ wouldn't be a reasonable assumption, hence $\mathbb E[Y_\mathrm{new}|x_\mathrm{new}]=x_\mathrm{new}\hat β$ wouldn't be the BLUP of $Y_\mathrm{new}$ and additionally, as you say, $\mathrm{Cov}(\varepsilon_\mathrm{new},\varepsilon)$ would have to be considered. Commented Aug 15, 2022 at 8:35