Goodmorning to List members,
I have a question about Regression Kriging.
Is it correct to define RK as interpolation method which combine regressive
interpolation + ordinary Kriging or regression residuals?
I compared GSIF package and gstat package but I have different
predictions...
I provide an example to explain better my question
-------------------------------------------------------------------------------------
library(gstat)
library(raster)
demo(meuse, echo=FALSE)
meuse
class(meuse)
names(meuse)
r<-raster(system.file("external/test.grd", package="raster"))
r<-as(r,'SpatialPixelsDataFrame') #to use predict function
names(r)<-'elev'
projection(r)<-projection(meuse)<-projection(meuse.grid)
plot(r,axes=T)
class(r)
summary(r)

####-------------------------Regression interpolation
fn<-lm(zinc~elev,meuse)
summary(fn)
gREG<-gstat(NULL,"zinc",zinc~elev,meuse)
pred.reg<-predict(gREG,r,progress='text')
summary(pred.reg)
plot(raster(pred.reg,'zinc.pred'))
plot(raster(pred.reg,'zinc.var'))

####-------------------------Universal Kriging interpolation
vu<-variogram(zinc~elev,meuse)
vu.fit<-fit.variogram(vu,vgm(1,"Exp",300,1))
plot(vu,vu.fit)
gUK<-gstat(NULL,"zinc",zinc~elev,meuse,model=vu.fit)
pred.uk<-predict(gUK,r,progress='text')
summary(pred.uk)
plot(raster(pred.uk,'zinc.pred'))
plot(raster(pred.uk,'zinc.var'))

####-------------------------Regression Kriging interpolation (1) {gstat}
vr<-variogram(residuals(fn)~1,meuse)
vr.fit<-fit.variogram(vr,vgm(1,"Exp",300,1))
plot(vr,vr.fit)
gRK<-gstat(NULL,"zinc",zinc~1,meuse,model=vr.fit)
pred.rk<-predict(gRK,r,progress='text')
pred.rk@data$zinc.pred<-pred.reg@data$zinc.pred+pred.rk@data$zinc.pred
pred.rk@data$zinc.var<-pred.reg@data$zinc.var+pred.rk@data$zinc.var #to be
analysed....
summary(pred.uk)
plot(raster(pred.uk,'zinc.pred'))
plot(raster(pred.uk,'zinc.var'))

####-------------------------Regression Kriging interpolation (2) {GSIF}
library(GSIF)
#?fit.gstatModel
Rv<-fit.gstatModel(meuse,zinc~elev,r,method='GLM',stepwise=T,vgmFun='Exp')
str(Rv)
summary(Rv@regModel)
Rv.fit<-Rv@vgmModel
class(Rv.fit)<-c("variogramModel","data.frame")
RKv<-variogram(zinc~1,meuse)
RKvu<-variogram(zinc~elev,meuse)
plot(x=RKvu$dist,y=RKvu$gamma,pch=20,xlab='distance',cex=1.1,ylab='gamma',ylim=c(0,max(RKv$gamma)))
points(x=RKv$dist, y=RKv$gamma,pch="+",cex=1.1,col="grey")
vline<-variogramLine(Rv.fit,maxdist=max(RKvu$dist),n=length(RKvu$dist))
lines(x=vline$dist,y=vline$gamma)
pred.rk0<-predict(Rv,r,nfold=nrow(meuse))
class(pred.rk0)
#plotKML(pred.rk0)
str(pred.rk0)
plot(raster(pred.rk0@predicted))
#
obj<-stack(raster(pred.uk),raster(pred.rk),raster(pred.rk0@predicted))
names(obj)<-c('UK_gstat','RK_gstat','RK_GSIF')
plot(obj)
pairs(obj) #correlation=1 but different histograms and not exactly the same
thing...
--------------------------------------------------------------------------------------------------------------------------------------

THANKS, maurizio

	[[alternative HTML version deleted]]

