6

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) 

# Load necessary packages
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)') 

classical variogram

# 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)') 

robust variogram

migrated from stats.stackexchange.com Aug 18 '14 at 14:09

This question came from our site for people interested in statistics, machine learning, data analysis, data mining, and data visualization.

  • 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? – Asad Ali Jan 29 '17 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... – adamdsmith Jan 30 '17 at 18:44
3

Interesting finding, and thanks for the perfectly prepared reproducible example! Both variograms you mention are meant to prepare for linear kriging in the next step, and you may not want that. The robust variogram was (IIRC) deviced for normal data with some pollution (outliers), but not for count data.

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).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.