Zev,
GLS only fits the trend; krige predicts trend + residual, unless you
specify BLUE=TRUE:
> theGLS$fit[1:5]
1 2 3 4 5
6.891951 6.703849 6.166895 5.873390 5.642714
> predict(g, meuse, BLUE=TRUE)[["var1.pred"]][1:5]
[generalized least squares trend estimation]
[1] 6.891951 6.703849 6.166895 5.873390 5.642714
Best regards,
On 10/28/2011 09:58 PM, Zev Ross wrote:
Hi All,
I was expecting the predictions from GSTAT to be the same as a GLS model
if the parameters and correlation structure are exactly the same but I'm
finding they're not. Why would that be? Below is an example using the
Meuse dataset.
Zev
data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y
theFormula <- "log(zinc)~sqrt(dist)"
lznr.vgm <- variogram(formula(theFormula), meuse)
lznr.fit <- fit.variogram(lznr.vgm, model = vgm(1, "Exp", 300, 1))
g <- gstat(formula=formula(theFormula), data = meuse, model = lznr.fit)
# get the coefficients so we can compare against GLS
lzn.coef<-predict(g, meuse, BLUE=c(TRUE,TRUE))
# generate the predictions using kriging with external drift
lzn.kriged <- predict(g, meuse, BLUE=FALSE)
# try to recreate what predict in GSTAT is doing using GLS
nug<-lznr.fit[1,2]
sill<-lznr.fit[1,2]+lznr.fit[2,2]
range<-lznr.fit[2,3]
# the parameters from theGLS should match exactly with the UKcoef
parameters
theGLS<-gls(formula(theFormula), data=meuse, na.action=na.omit,
correlation=corExp(c(range, nug/sill),
form=~x+y, nugget=TRUE,fixed=TRUE), method="ML")
# confirm that the coefficients are the same
lzn.coef
summary(theGLS)$coef
# yes they are identical
> lzn.coef
coordinates var1.pred var1.var
(Intercept) (181072, 333611) 6.985991 0.02507861
sqrt(dist) (181072, 333611) -2.551849 0.07302336
> summary(theGLS)$coef
(Intercept) sqrt(dist)
6.985991 -2.551849
# the predictions/fits are different
> log(meuse@data$zinc[1:5]) #true data
[1] 6.929517 7.039660 6.461468 5.549076 5.594711
> lzn.kriged@data$var1.pred[1:5] # output from kriging model
[1] 6.929517 7.039660 6.461468 5.549076 5.594711
> theGLS$fit[1:5] # output from GLS
6.891951 6.703849 6.166895 5.873390 5.642714
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics e.pebe...@wwu.de
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
