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?
 A: "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).
A: 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.
A: The other answers are good I just wanted to add something about 'accounting for' spatial autocorrelation.  Sometimes this claim is made more strongly along the lines of "accounting for spatial autocorrelation not explained by the covariates".
This can present a misleading picture of what the spatial smooth does.  It is not like there is some orderly queue in the likelihood where the smooth patiently waits for the covariates to go first and then smooth will mop up the 'unexplained' parts.  In reality they all get a chance to explain the data.  
This paper with an aptly named title presents the issue really clearly, although it is from the point of view of a CAR model the principles apply to GAM smooths.
Adding Spatially-Correlated Errors Can Mess Up the Fixed Effect You Love
The 'solution' in the paper is to smooth the residuals instead of smoothing on space.  That would have the effect of allowing your covariates to explain what they can. Of course, there are many applications in which this would not be a desireable solution. 
A: Spatial correlation is simply how the x and y coordinates relate to the magnitude of the resulting surface in the space.  So the autocorrelation between the coordinates can be expressed in terms of a functional relationship between the neighboring points.
