# Classical or robust variogram for incorporation into generalized linear model

I'm modeling counts of organisms over time at eleven locations. I'd like to account for temporal autocorrelation in the counts, assuming it's present.

As the data are not equally spaced in time, I'm exploring empirical variograms (and comparing variogram model fits) based on the residuals of a generalized linear model assuming independence.

As shown below, the variogram based on the classical semivariance estimator looks quite different from the robust (Cressie) estimator. Furthermore, a fit of most variogram models to the classical variogram produces singularity errors (with good reason based on the variogram plot) while a decent exponential model (e.g.) can be fitted to the robust variogram. Unfortunately, it seems as though `nlme:::lme` (which I'm using for the generalized model) only uses the classical variogram when fitting correlation structures.

Thus, my questions:

(1) Is the use of the robust estimator for constructing the correlation structure justified in this case?
(2) If so, is it appropriate to fit the the robust variogram model in, e.g., `gstat` or `geoR`, and then specify/fix the correlation structure in the `lme` fit?

``````# Download data
datURL <- "https://www.dropbox.com/s/54q9ocwztt3swap/rsiggeo_variogram.csv"
dat <- repmis::source_data(datURL, sep = ",", header = TRUE)

library(gstat)
library(sp)

# Convert to appropriate class
# Y coordinates are arbitrarily assigned values for each of 11 sites
# (count locations; assumed independent); Y coords were assigned such
# that each site is separated from others by more than the cutoff distance
# used in variogram construction
coordinates(dat) <- ~x+y

# Classic variogram model with example exponential fit
varClassic <- variogram(resid ~ 1, dat, cutoff = 180, width = 180/25)
# Note singularity error; starting values determine 'fit'
expClassic <- fit.variogram(varClassic, model = vgm(1, "Exp", 1, 0))
plot(varClassic, expClassic, xlab = 'Distance (days)')
``````

``````# Robust (Cressie) variogram model with example exponential fit
varCressie <- variogram(resid ~ 1, dat, cutoff = 180, width = 180/25, cressie=TRUE)
expCressie <- fit.variogram(varCressie, model = vgm(0.3, "Exp", 50, 0.1))
plot(varCressie, expCressie, xlab = 'Distance (days)')
``````

• I had a look of your data and I am confused about your coordinates x y. When you take time as coordinate then does not it become one dimensional interval? I have been doing this by taking my x-coordinate as zero and y-coordinates as the time (hour, day, week or month etc). In this case, when you compute the euclidean distance for a point pair, it is just the time lap between those two points. Thus the distance becomes a time lap. Or, if I am mistaking it and doing it wrong, please guide me how do I deal with time? Commented Jan 29, 2017 at 13:43
• You could do it that way (i.e., by assigning all the y coordinates to zero), but then the distinction among sites is lost and the variogram would be estimated from the entirety of the data. In this case, that would distinctly muddle the situation because sites differ considerably in counts... Commented Jan 30, 2017 at 18:44

I would advice to look at model-based geostatistics, and since you mentioned count data, you'd probably end up using package `geoRglm`. The canonical reference is the book with the same title by Ribeiro and Diggle. You should also look at package `geoCount` which was just published in issue 11 found here; the paper includes comparisons with `geoRglm` and `INLA` (yet another alternative).