[R-sig-Geo] How to account for spatial autocorrelation
Luciano La Sala
lucianolasala at yahoo.com.ar
Wed Apr 17 16:29:34 CEST 2013
Dear experts,
I have a dataset which consists of 450 observations (cattle rising ranches)
from a survey made in two different provinces. Each ranch was classified as
diseases (1) / not diseased (0) based on diagnostic testing for certain
health event made only on some fraction of the cattle from each ranch.
My independent variables are 1) total heads in ranch, 2) total cattle in
ranch. My data looks like this:
Department Province Lat Long Total_heads Total_cattle
Disease_status
Atreuco La Pampa -36.97 -63.73 1900 600 0
Atreuco La Pampa -37.19 -63.69 2340 200 0
Atreuco La Pampa -37.20 -63.45 3200 150 1
BuenosAires BuenosAires -37.18 -63.44 1000 230 1
BuenosAires BuenosAires -37.08 -64.03 2300 85 0
Starting from a fixed-effects, only-intercept logistic model, none of the
independent variables improved the model fit, which leaves me with a null
model as the best fitting one.
Looking at a semivariogram of the standardized residuals versus distance,
model residuals seem to be autocorrelated. A missing covariate may the
reason behind residual autocorrelation. To prove the existence of residual
autocorrelation, I want to add different covariance structures to the null
model and compare the likelihoods of these models. As far as I know, adding
spatial correlation structures requires using the nlme package, and
the lme command needs a grouping variable (which my data lack). Therefore, I
created a dummy variable that has the same value for all 450 observations
and specified my mixed-effect null model as follows:
dummy <- rep(1, 450)
my.data <- cbind(my.data, dummy)
null.model <- lme(fixed = Disease_status ~ 1, data = my.data, random = ~ 1 |
dummy, method = "ML")
summary(null.model)
Linear mixed-effects model fit by maximum likelihood
Data: Data.new
AIC BIC logLik
454.3587 466.6798 -224.1793
Random effects:
Formula: ~1 | dummy
(Intercept) Residual
StdDev: 6.979292e-06 0.3986575
Fixed effects: Estatus_campo ~ 1
Value Std.Error DF t-value p-value
(Intercept) 0.1982183 0.0188348 448 10.52405 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.4972145 -0.4972145 -0.4972145 -0.4972145 2.0112046
Number of Observations: 449
Number of Groups: 1
The above model without a covariance structure would represent the baseline
likelihood, which I hoped to improve with information about the spatial
structure of the data. When I tried to update my null model to include an
exponential spatial correlation structure and see if the model with spatial
correlation fits better than one without, an error message pops up:
exp.sp <- update(null.model, correlation = corExp(1, form = ~ Long + Lat),
method = "ML")
summary(exp.sp)
Error in getCovariate.corSpatial(object, data = data) : cannot have zero
distances in "corSpatial"
Question: what is the cause behind this error and how can I solve it?
Thank you very much in advance!
Luciano La Sala
More information about the R-sig-Geo
mailing list