[R-sig-Geo] GSTAT predictions vs GLS

[Edzer Pebesma]

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 at data$zinc[1:5]) #true data
> [1] 6.929517 7.039660 6.461468 5.549076 5.594711
>
>  > lzn.kriged at 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.pebesma at wwu.de