Dear Lyndon, Have a look at the normalised residuals when creating a variogram. Those take the correlation structure into account. And have a look at the parameters of the correlation structure. Are they from a similar magnitude as you would expect from the OLS variogram? If not, rerun the GLS with starting values for range and sill based on the OLS variogram. I had some strange results in the past with range parameters being smaller than the smallest distance between two points. Supplying reasonable starting values yielded better results.
HTH, Thierry ------------------------------------------------------------------------ ---- ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek team Biometrie & Kwaliteitszorg Gaverstraat 4 9500 Geraardsbergen Belgium Research Institute for Nature and Forest team Biometrics & Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey > -----Oorspronkelijk bericht----- > Van: r-sig-geo-boun...@stat.math.ethz.ch > [mailto:r-sig-geo-boun...@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. > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
