-----Oorspronkelijk bericht-----
Van: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] Namens Lyndon Estes
Verzonden: maandag 14 juni 2010 18:15
Aan: r-sig-geo
Onderwerp: [R-sig-Geo] Advice on GLS on residual variogram?
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.