[R-sig-Geo] GSTAT predictions vs GLS

[Zev Ross]

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



-- 
Zev Ross
ZevRoss Spatial Analysis
120 N Aurora, Suite 3A
Ithaca, NY 14850
607-277-0004 (phone)
866-877-3690 (fax, toll-free)
zev at zevross.com