[R-sig-Geo] Advice on GLS on residual variogram?

[Lyndon Estes]

Hello,

I am attempting to implement regression kriging, following methods
recommended by Hengl et al. (2007; and others), for % soil organic
carbon over South Africa. I was hoping to ask for advice as to whether
the results I am getting from GLS make sense.

For background, I have 3300+ soil profiles providing A horizon OC %,
and I have derived 8 spatial predictors including slope, solar
radiation, a topographic moisture index, etc. These have been
transformed using principal components analysis (in ArcGIS).

My question concerns variograms resulting from the first component of
the methodology, which is to find the coefficients and residuals of a
GLS model as follows:

1. Use OLS to fit my model:

oc<-read.dbf('~/oc/ocdat.dbf')
oc_2<-as.data.frame(oc[,c(2,5:6,21:27,207:225)])

gls.all<-gls(log(CTOP)~PCB1+PCB2+PCB3+PCB4+PCB5+PCB6+PCB8,data=oc_2)
# Note: GLS is OLS if correlation structure is not specified

2. Find the GLS coefficients (and residuals) using an appropriate
autocorrelation structure. In this case, fitting a variogram to the
OLS residuals suggested a spherical autocorrelation structure:

gls.all.update<-update(gls.all,correlation=corSpher(form=~X+Y,nugget=TRUE))

Step 2 took a long time to complete, given my dataset--an overnight
run (not sure how many hours though) using 64-bit R2.10.1 on a MacBook
Pro with a 2.4 GHZ processor and 4 GB of RAM. I am not sure that the
results make sense, however, as the GLS shows greater autocorrelation
in the residuals then the original OLS residuals. The following
produced the illustration of the plotted residual variograms posted
here (http://sites.google.com/site/ldemisc/variogram):

# Fit variograms to residuals of OLS and GLS
oc.var<-variogram(residuals(gls.all)~1,oc_2)
oc.var.update<-variogram(residuals(gls.all.update)~1,oc_2)

# Create variograms plots
oc.var.update.pl<-plot(oc.var.update,main="GLS residuals")
oc.var.pl<-plot(oc.var,main="OLS residuals")

# Display variograms side-by-side
oc.var.pl$x.limits<-oc.var.update.pl$x.limits
oc.var.pl$y.limits<-oc.var.update.pl$y.limits
print(oc.var.pl,split=c(1,1,2,1),more=TRUE)
print(oc.var.update.pl,split=c(2,1,2,1),more=FALSE)

Does it seem sensible that GLS residuals show a stronger degree of
spatial autocorrelation (with no sign of a sill) than OLS? For
comparison with another spatially autocorrelated dataset, I used the
meuse dataset (following an example from Hengl 2009) with the models:

data(meuse)
coordinates(meuse)=~x+y
meu.ols<-gls(log(zinc)~dist.m+ffreq,meuse)
meu.gls<-update(meu.ols,correlation=corExp(form=~x+y))

Plotting the variograms:

zinc.vgmOLS<-variogram(residuals(meu.ols)~1,meuse)
ols.vgm.pl<-plot(zinc.vgmOLS,main="OLS plot")

zinc.vgm.GLS<-variogram(residuals(meu.gls)~1,meuse)
gls.vgm.pl<-plot(zinc.vgm.GLS,main="GLS plot")

# Display variograms side-by-side
ols.vgm.pl$x.limits<-gls.vgm.pl$x.limits
ols.vgm.pl$y.limits<-gls.vgm.pl$y.limits
print(ols.vgm.pl,split=c(1,1,2,1),more=TRUE)
print(gls.vgm.pl,split=c(2,1,2,1),more=FALSE)

In contrast, this shows nearly identical spatial autocorrelation
patterns for OLS and GLS.

I would greatly appreciate any advice regarding the (seemingly)
unusual variogram results, and/or clearing up of any misunderstandings
I might have.

Thanks in advance.

Regards, Lyndon


Hengl, T., G.B.M. Heuvelink, and D.G. Rossiter. 2007. About
regression-kriging: From equations to case studies. Computers &
Geosciences 33, no. 10 (October): 1301-1315.

Hengl, T. 2009. A Practical Guide to Geostatistical Mapping.
Luxembourg: Office for Official Publications of the European
Communities.