[R-sig-Geo] How to account for spatial autocorrelation

[Luciano La Sala]

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