I'm fitting a non-stationary model using mgcv (family: gaulss()) where the data have been collected at different points in space. The dataset is large (N>10000) and was collected at essentially equal distances from each other across the spatial domain. Here's a simplified version of what the model looks like in R:
f <- y ~ s(x) #Mean model
f2 <- ~ s(x) #Variance model
flist <- list(f, f2) #List of model formulae
mod <- gam(flist, family=gaulss(), data=dat) #Fit location-scale (non-stationary) model
Simple residual diagnostics (Q-Q plots) appear OK, but since the data were collected over space, I'd like to look at spatial autocorrelation in the residuals. Using the variogram function from the gstat library, I've made variograms of two types of residuals: response residuals ($y - \hat{y}$) and "deviance" residuals ($\frac{y-\hat{y}}{\sigma}$). As you can see, the results are very different:
My question is: which type of residual is appropriate for assessing spatial (or temporal) autocorrelation in non-stationary models?
