Formula for prediction interval for spatial regression

Question

I have a spatial error model that I have estimated and would like to create a figure with the prediction interval. By spatial error model, I mean that the errors are spatially correlated and thus OLS is unbiased but inefficient (not a spatially auto-regressive model).

I used the geostatistical approach to spatial statistics and estimated a semivariogram for my data and used the semivariogram estimates to create a weighting matrix to account for the spatial autocorrelation in the data. I think the answer to this question though would help anyone using weights in their regression regardless of the exact structure of their weighting matrix.

I see from this question that the formula for a 95% prediction interval for OLS is
$$ \hat{y} \pm 1.96 \hat{\sigma} \sqrt{1 + \mathbf{X}^* (\mathbf{X}'\mathbf{X})^{-1} (\mathbf{X}^*)'}. $$

Given my limited intuition, my guess would be that the prediction interval for the case with weights would be this (where $\Omega$ is the weighting matrix):

$$ \hat{y} \pm 1.96 \hat{\sigma} \sqrt{1 + \mathbf{X}^* (\mathbf{X}'\mathbf{\Omega}^{-1}\mathbf{X})^{-1} (\mathbf{X}^*)'} $$

and was hoping that I could get some feedback from the community here about the correctness or if I am even in the right ballpark.

Answer 1

Note that @whuber adresses all the points raised in the question, I will provide only the mathematical details. The regression with spatial error assumes that model is the following:

$$y=X\beta+u$$

where $Euu'=\Omega$. The efficient estimate of $\beta$ then is

$$\hat\beta=(X'\Omega X)^{-1}X'\Omega y$$

When you have new observation it also follows the model

$$y^*=X^*\beta+u^*$$

We predict $y^*$ by $X^*\beta$ so the prediction error is

$$y^*-\hat{y^*}=X^*\hat\beta+u^*$$

For linear model $u^*$ will be independent of $u$ so

$$Var(y^*-\hat{y^*})=Var(X^*\hat\beta+u^*)=Var(X^*\hat\beta)+Var(u^*)=\sigma^2+\sigma^2X^{* } (X'X)^{-1})(X^{* })'$$

For the spatial error model we have \begin{align*} Var(X^{* }\hat\beta)= Var(X^{* } (X'\Omega X)^{-1}X'\Omega u) =X^{* } (X'\Omega X)^{-1})(X^{* })' \end{align*}

But now $u^*$ correlates with $u$ and we do not have such convenient decomposition. Thus the formula in your question cannot be used.

The general approach how to calculate the error is describe in this article. It shows how to get the best unbiased linear prediction in linear model when the model errors are correlated. Note that this heavily relies on the fact that you have reliable estimate of $\Omega$, so @whuber answer applies very strongly. Even if you have correct formulas they might give bad results if it is known that $\Omega$ cannot be reliably estimated.

Answer 2

Assuming $\Omega$ has been estimated with no appreciable error--which is rarely the case--the correct formula is given by the equation for the kriging prediction error. It needs to replace everything under the square root sign. What is missing from the second formula you give is any explicit accounting of the covariance between the value at the estimation location and the data values.

Some recent approaches, such as presented in Diggle & Ribeiro's "Model-based Geostatistics" (implemented as an R package) have the additional merit of incorporating the estimation error of $\Omega$ in the prediction uncertainty.