# Why does including latitude and longitude in a GAM account for spatial autocorrelation?

I have produced generalized additive models for deforestation. To account for spatial-autocorrelation, I have included latitude and longitude as a smoothed, interaction term (i.e. s(x,y)).

I've based this on reading many papers where the authors say 'to account for spatial autocorrelation, coordinates of points were included as smoothed terms' but these have never explained why this actually accounts for it. It's quite frustrating. I've read all the books I can find on GAMs in the hope of finding an answer, but most (e.g. Generalized Additive Models, an Introduction with R, S.N. Wood) just touch on the subject without explaining.

I'd really appreciate it if someone could explain WHY the inclusion of latitude and longitude accounts for spatial autocorrelation, and what 'accounting' for it really means - is it simply enough to include it in the model, or should you compare a model with s(x,y) in and a model without? And does the deviance explained by the term indicate the extent of spatial autocorrelation?

-
If it is relevant, I used the 'bam' function from the 'mgcv' package in R. –  gisol Sep 1 '12 at 14:05
Also, I've tested for spatial autocorrelation using Moran's I. –  gisol Sep 1 '12 at 14:19
–  Macro Sep 1 '12 at 15:29
Given the answers here, we might flag the other Q @Macro links to as a duplicate of this one so people coming across that one see the Answers here, especially that of whuber. –  Gavin Simpson Sep 1 '12 at 17:20
+1 @GavinSimpson - by the way, note that you do have the power to cast close votes, enough of which will lead to the two questions being merged. –  Macro Sep 1 '12 at 20:11
show 1 more comment

The main issue in any statistical model is the assumptions that underlay any inference procedure. In the sort of model you describe, the residuals are assumed independent. If they have some spatial dependence and this is not modelled in the sytematic part of the model, the residuals from that model will also exhibit spatial dependence, or in other words they will be spatially autocorrelated. Such dependence would invalidate the theory that produces p-values from test statistics in the GAM for example; you can't trust the p-values because they were computed assuming independence.

You have two main options for handling such data; i) model the spatial dependence in the systematic part of the model, or ii) relax the assumption of independence and estimate the correlation between residuals.

i) is what is being attempted by including a smooth of the spatial locations in the model. ii) requires estimation of the correlation matrix of the residuals often during model fitting using a procedure like generalised least squares. How well either of these approaches deal with the spatial dependence will depend upon the nature & complexity of the spatial dependence and how easily it can be modelled.

In summary, if you can model the spatial dependence between observations then the residuals are more likely to be independent random variables and therefore not violate the assumptions of any inferential procedure.

-
Thanks for your clear answer Gavin. What makes spatial autocorrelation fundamentally different from any gradient not included in the model? Say your study area was on a sloping hill, and the species of interest preferred lower habitat to higher habitat. Failing to include the elevation in the model would leave a structure in the residuals, would it not? Is it simply that spatial autocorrelation is (or was) forgotten about or not considered? (P.S. perhaps this is a poor example as inclusion of lat,long would account for this effect too). –  gisol Sep 1 '12 at 16:43
Yes. I suspect that in the examples you have looked at either the spatial component was of interest so was modelled explicitly via a smooth of lat/lon or the spatial component was a nuisance term but needed to be modelled to leave the residuals i.i.d. If the "spatial" component is better modelled via a different variable (e.g. elevation in you comment) then a smooth of that variable would be used instead of the spatial locations. –  Gavin Simpson Sep 1 '12 at 17:03

"Spatial autocorrelation" means various things to various people. An overarching concept, though, is that a phenomenon observed at locations $\mathbf{z}$ may depend in some definite way on (a) covariates, (b) location, and (c) its values at nearby locations. (Where the technical definitions vary lie in the kind of data being considered, what "definite way" is postulated, and what "nearby" means: all of these have to be made quantitative in order to proceed.)

To see what might be going on, let's consider a simple example of such a spatial model to describe the topography of a region. Let the measured elevation at a point $\mathbf{z}$ be $y(\mathbf{z})$. One possible model is that $y$ depends in some definite mathematical way on the coordinates of $\mathbf{z}$, which I will write $(z_1,z_2)$ in this two-dimensional situation. Letting $\varepsilon$ represent (hypothetically independent) deviations between the observations and the model (which as usual are assumed to have zero expectation), we may write

$$y(\mathbf{z}) = \beta_0 + \beta_1 z_1 + \beta_2 z_2 + \varepsilon(\mathbf{z})$$

for a linear trend model. The linear trend (represented by the $\beta_1$ and $\beta_2$ coefficients) is one way to capture the idea that nearby values $y(\mathbf{z})$ and $y(\mathbf{z}')$, for $\mathbf{z}$ close to $\mathbf{z}'$, should tend to be close to one another. We can even calculate this by considering the expected value of the size of the difference between $y(\mathbf{z})$ and $y(\mathbf{z}')$, $E[|y(\mathbf{z}) - y(\mathbf{z}')|]$. It turns out the mathematics is much simpler if we use a slightly different measure of difference: instead, we compute the expected squared difference:

\eqalign{ E[\left(y(\mathbf{z}) - y(\mathbf{z}')\right)^2] &= E[\left(\beta_0 + \beta_1 z_1 + \beta_2 z_2 + \varepsilon(\mathbf{z}) - \left(\beta_0 + \beta_1 z_1' + \beta_2 z_2' + \varepsilon(\mathbf{z}')\right)\right)^2] \\ &=E[\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)' + \varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)^2] \\ &=E[\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)^2 \\ &\quad+ 2\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)\left(\varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)\\ &\quad+ \left(\varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)^2] \\ &=\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)^2 + E[\left(\varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)^2] }

This model is free of any explicit spatial autocorrelation, because there is no term in it directly relating $y(\mathbf{z})$ to nearby values $y(\mathbf{z}')$.

An alternative, different, model ignores the linear trend and supposes only that there is autocorrelation. One way to do that is through the structure of the deviations $\varepsilon(\mathbf{z})$. We might posit that

$$y(\mathbf{z}) = \beta_0 + \varepsilon(\mathbf{z})$$

and, to account for our anticipation of correlation, we will assume some kind of "covariance structure" for the $\varepsilon$. For this to be spatially meaningful, we will assume the covariance between $\varepsilon(\mathbf{z})$ and $\varepsilon(\mathbf{z}')$, equal to $E[\varepsilon(\mathbf{z})\varepsilon(\mathbf{z}')]$ because the $\varepsilon$ have zero means, tends to decrease as $\mathbf{z}$ and $\mathbf{z}'$ become more and more distant. Because the details do not matter, let's just call this covariance $C(\mathbf{z}, \mathbf{z}')$. This is spatial autocorrelation. Indeed, the (usual Pearson) correlation between $y(\mathbf{z})$ and $y(\mathbf{z}')$ is

$$\rho(y(\mathbf{z}), y(\mathbf{z}')) = \frac{C(\mathbf{z}, \mathbf{z}')}{\sqrt{C(\mathbf{z}, \mathbf{z})C(\mathbf{z}', \mathbf{z}')}}.$$

In this notation, the previous expected squared difference of $y$'s for the first model is

\eqalign{ E[\left(y(\mathbf{z}) - y(\mathbf{z}')\right)^2] &= \left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)^2 + E[\left(\varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)^2] \\ &=\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)^2 + C_1(\mathbf{z}, \mathbf{z}) + C_1(\mathbf{z}', \mathbf{z}') }

(assuming $\mathbf{z} \ne \mathbf{z}'$) because the $\varepsilon$ at different locations have been assumed to be independent. I have written $C_1$ instead of $C$ to indicate this is the covariance function for the first model.

When the covariances of the $\varepsilon$ do not vary dramatically from one location to another (indeed, they are usually assumed to be constant), this equation shows that the expected squared difference in $y$'s increases quadratically with the separation between $\mathbf{z}$ and $\mathbf{z}'$. The actual amount of increase is determined by the trend coefficients $\beta_0$ and $\beta_1$.

Let's see what the expected squared differences in the $y$'s is for the new model, model 2:

\eqalign{ E[\left(y(\mathbf{z}) - y(\mathbf{z}')\right)^2] &= E[\left(\beta_0 + \varepsilon(\mathbf{z}) - \left(\beta_0 + \varepsilon(\mathbf{z}')\right)\right)^2] \\ &=E[\left(\varepsilon(\mathbf{z}) - \varepsilon(\mathbf{z}')\right)^2] \\ &=E[\varepsilon(\mathbf{z})^2 - 2 \varepsilon(\mathbf{z})\varepsilon(\mathbf{z}') + \varepsilon(\mathbf{z}')^2] \\ &=C_2(\mathbf{z}, \mathbf{z}) - 2C_2(\mathbf{z}, \mathbf{z}') + C_2(\mathbf{z}', \mathbf{z}'). }

Again this behaves in the right way: because we figured $C_2(\mathbf{z}, \mathbf{z}')$ should decrease as $\mathbf{z}$ and $\mathbf{z}'$ become more separated, the expected squared difference in $y$'s indeed goes up with increasing separation of the locations.

Comparing the two expressions for $E[\left(y(\mathbf{z}) - y(\mathbf{z}')\right)^2]$ in the two models shows us that $\left(\beta_1 (z_1-z_1') + \beta_2 (z_2-z_2)'\right)^2$ in the first model is playing a role mathematically identical to $-2C_2(\mathbf{z}, \mathbf{z}')$ in the second model. (There's an additive constant lurking there, buried in the different meanings of the $C_i(\mathbf{z}, \mathbf{z})$, but it doesn't matter in this analysis.) Ergo, depending on the model, spatial correlation is typically represented as some combination of a trend and a stipulated correlation structure on random errors.

We now have, I hope, a clear answer to the question: one can represent the idea behind Tobler's Law of Geography ("everything is related to everything else, but nearer things are more related") in different ways. In some models, Tobler's Law is adequately represented by including trends (or "drift" terms) that are functions of spatial coordinates like longitude and latitude. In others, Tobler's Law is captured by means of a nontrivial covariance structure among additive random terms (the $\varepsilon$). In practice, models incorporate both methods. Which one you choose depends on what you want to accomplish with the model and on your view of how spatial autocorrelation arises--whether it is implied by underlying trends or reflects variations you wish to consider random. Neither one is always right and, in any given problem, it's often possible to use both kinds of models to analyze the data, understand the phenomenon, and predict its values at other locations (interpolation).

-
+1 - it's nice to see the link between two approaches for handling spatial dependence. Great answer, whuber! –  Macro Sep 1 '12 at 16:59
Very comprehensive, thank you. It will take me a few moments to think all this through. –  gisol Sep 1 '12 at 17:04
+1 Simply marvellous! –  Gavin Simpson Sep 1 '12 at 17:14
If all statistical writing were of this ilk there would be a lot more clear-thinking applied statistical work in the world. Beautifully done. –  Ari B. Friedman Mar 11 at 12:21
+1 wow, another great answer by whuber. Very helpful! –  TaylerJones Dec 3 at 17:28