[R-sig-ME] modeling spatial autocorrelation and heteroscedascity with gls / nlme

[Guido Lorenz]

Dear R users,

my question is about groupwise modeling of spatial autocorrelation in mixed effect models: In a soil study, five different sites are to be compared respect to specific soil properties. A regular sampling grid with 32 points in each site was applied. Due to heteroscedascity and a supposed spatial correlation of the data, a mixed effects model was formulated, with the "site" (ID.sitio) as a fixed effect, and different variances and spatial autocorrelation, also stratified by the same factor, "site", using the gls module of nlme package: 

> lmm.dap3a <- gls(Dap.sa ~ 1 + ID.sitio, Pd2006mm, correlation=corSpher(form=~(easting.m+northing.m)|ID.sitio, nugget=TRUE, metric="euclidean"), weights=varIdent(form=~1|ID.sitio), na.action=na.omit, method="REML") 

>  summary(lmm.dap3a)
Generalized least squares fit by REML
  Model: Dap.sa ~ 1 + ID.sitio 
  Data: Pd2006mm 
        AIC       BIC   logLik
  -280.4552 -243.9341 152.2276

Correlation Structure: Spherical spatial correlation
 Formula: ~(easting.m + northing.m) | ID.sitio 
 Parameter estimate(s):
     range     nugget 
17.7285895  0.2560010 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | ID.sitio 
 Parameter estimates:
     pag1      pag2      pag3      pag4      pag5 
1.0000000 1.4464638 0.6195818 1.0094856 0.7684725 

Coefficients:
                  Value  Std.Error  t-value p-value
(Intercept)   1.3930856 0.03141931 44.33851  0.0000
ID.sitiopag2  0.0520236 0.05601737  0.92871  0.3545
ID.sitiopag3 -0.1816564 0.03708524 -4.89835  0.0000
ID.sitiopag4  0.0160424 0.04424861  0.36255  0.7174
ID.sitiopag5 -0.1334217 0.03983743 -3.34915  0.0010

 Correlation: 
             (Intr) ID.st2 ID.st3 ID.st4
ID.sitiopag2 -0.561                     
ID.sitiopag3 -0.847  0.475              
ID.sitiopag4 -0.710  0.398  0.602       
ID.sitiopag5 -0.789  0.442  0.668  0.560

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-2.24247234 -0.61313647 -0.01625797  0.67134893  2.21471264 

Residual standard error: 0.1117062 
Degrees of freedom: 160 total; 155 residual



Nevertheless, the output shows that the variances were described as group-specific, but the spatial correlation is expressed as a unique semivariogram for all the sites. Supposing that spatial correlation may be also specific for each sampling unit, how can one specify the model to account for this, using the gls or or some other algorithm?

Thanks in advance for any advice,

Guido Lorenz