Modeling a spatial trend by regression with the $(x,y)$ coordinates as predictors

Question

I plan to include coordinates as covariates in the regression equation in order to adjust for the spatial trend that exists in the data. After that, I want to test residuals on spatial autocorrelation in random variation. I have several questions:

1. Should I perform linear regression in which only independent variables are $x$ and $y$ coordinates and then tests residuals on spatial autocorrelation, or should I rather include not only coordinates as covariates but also other variables and then test residuals.

2. If I expect to have quadratic trend, and then include not only $x,y$, but also , $xy$, $x^2$ and $y^2$, but then some of them ($xy$ and $y^2$) have the $p$-value higher than the threshold --should I exclude those variables with higher $p$-value as being nonsignificant? How should I then interpret the trend, it is certainly not quadratic anymore?

3. I guess I should treat $x$ and $y$ coordinates as any other covariates, and test them on having linear relationship with dependent variable by constructing partial residual plots ... but then once I transform them (if they show they need transformation), that will not be that kind of trend any more (especially if I include $xy$, $x^2$ and $y^2$ for quadratic trend). It may show that $x^2$, for example, needs transformation, while $x$ does not or so? How should I react in these situations?

Thank you.

Answer 1

I think you might be better off fitting a linear mixed effects model with spatially-correlated random effects (sometimes called a geostatistical model). Assuming your data is Gaussian, you specify a model of the form:

$ Y_i = \mu_i + S_i + \epsilon_i, $

for $n$ observations $1 \leq i \leq n$, with $\epsilon \sim N(0,\tau^2)$ representing iid errors and $\mathbb{S} \sim MVN(\mathbb{0},\sigma^2 R)$ representing your spatial terms (where $\mathbb{S} = \{S_1,...,S_n\}$). The mean $\mu_i$ could be a function of other covariates (i.e. $\mu_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}$ etc.) or it could be just be a constant (might be best to start with the latter for simplicity).

The correlation matrix $R$ for the spatial terms (which determines how correlated you think each observation should be) can be specified by looking at the empirical variogram. Generally the correlation between observations is chosen to only depend on the distance between them (this is where your coordinates come into the model).

Chapter 2 of Model-based geostatistics by Diggle and Ribeiro (2000) should give you a more detailed introduction. The R package geoR has many procedures for fitting geostatistical models, so you may find it useful (see http://cran.r-project.org/web/packages/geoR/geoR.pdf).