[R-sig-Geo] Problem in regression kriggin

[Ashton Shortridge]

Dear Abdul,

there are a couple of strange things about your code (the cross-variogram 
section actually performs cross-validation), but there are two problems with 
your code for the regression kriging in particular:
1. in my version of R at least, specifying 'data=data' throws the error you 
describe. This way seems to work better:
prk.rich<-krige(Richness~Elev, data, newdata=grd,model=fit.sph1)

2. More critically, the prediction object you wish to predict to (grd) must 
have the independent variable(s) in the model as attributes; in your case, 
there is one variable called Elev. However, you do not have that. If you could 
import a DEM and use that as the prediction object, this might work. For 
example, if I change your code to ordinary kriging, it works fine:
prk.rich<-krige(Richness~1, data, newdata=grd,model=fit.sph1)

though of course this is the wrong variogram model to use in ordinary kriging!

For more on regression kriging, check out the following wiki page, especially 
the R+gstat section:
http://spatial-analyst.net/wiki/index.php?title=Regression-kriging_guide

On 01/08/13, halim10-fes, wrote:
> Dear Group,
> 
> Happy New Year to Everybody.
> 
> I am a novice user of gstat and geostatistics. I have a data set (please
> see attached: data.csv) containing coordinates (East_x,
> North_y),Elev(Elevation), and Richness (Plant Species) collected from 125
> plots of 10m x 10m plots. I want to interpolate plant species richness
> from the relationship between Richness and Elevation of the study area
> using Regression Krigging. I have written the following codes:
> 
> #Variogram of residuals
> 
#==========================================================================
> =
> 
> data<-read.csv("E://data.csv",header=TRUE)
> attach(data)
> names(data)
> m2=glm(Richness~Elev)
> summary(m2)
> residuals<-m2$residuals
> data$Residuals<-residuals
> 
> library(maptools)
> 
> coordinates(data) = ~ East_x+North_y
> 
> library(gstat)
> 
> var <- variogram(Residuals~1,
> data=data,alpha=c(0,22.5,45,67.5,90,112.5,135),
> tol.hor=30,cutoff=1000,width=67)
> plot(var)
> 
> #Variogram shows that major axis is 45 & minor 112.5; Fitting variogram in
> #major axis
> 
> var.45<-var$dir.hor==45
> plot(var[var.45,])
> 
> eye.sph<-vgm(psill=10.5,model="Sph",range=705,nugget=2)
> plot(var[var.45,],eye.sph)
> 
> 
> #Fit anisotropic variogram
> 
> fit.sph<-fit.variogram(var,model=vgm(psill=10.5,model="Sph",range=705,
>                                      nugget=2,anis=c(45,0.72)))
> plot(var, model=fit.sph, as.table=TRUE)
> 
> fit.exp<-fit.variogram(var, model=vgm(psill=10.5,model='Exp',range=705/3,
>                                       nugget=2, anis=c(45, 0.72)))
> plot(var, model=fit.exp, as.table=TRUE)
> 
> fit.gau<-
fit.variogram(var,model=vgm(psill=10.5,model='Exp',range=705/sqrt(
> 3), nugget=2, anis=c(45, 0.72)))
> plot(var, model=fit.gau, as.table=TRUE)
> 
> #Cross Variogram
> 
> cv.sph<-krige.cv(Residuals~1,data,fit.sph)
> View(cv.sph)
> mean(cv.sph$zscore) #-0.0001968425
> var(cv.sph$zscore)  #1.114187
> 
> cv.exp<-krige.cv(Residuals~1,data,fit.exp)
> View(cv.exp)
> mean(cv.exp$zscore)  #-1.243423e-05
> var(cv.exp$zscore)   #1.116687
> 
> cv.gau<-krige.cv(Residuals~1,data,fit.gau)
> View(cv.gau)
> mean(cv.gau$zscore)  #-1.235693e-05
> var(cv.gau$zscore)   #1.116684
> 
> 
> #Cross variogram shows that spherical model shows the best fit. The sill
> varies in different direction of the varigram. So, we will have to fit a
> zonal #anisotropic variogram with spherical model (fit.sph). Variogram
> shows that #22.5 deg has the highest sill and 112.5 deg should have the
> lowest sill. We #will fit the zonal variogram to the 112.5 deg direction
> with an anisotropic #ratio of 293/(1600*10000)
> 
> #Fitting zonal anisotropic variogram
> 
> fit.sph1<-vgm(psill=2.5,"Sph",range=1600*10000,anis=c(122.5,1.825e-05),
>               add.to=fit.sph)
> plot(var,fit.sph1,as.table=TRUE)
> 
> #Cross Variogram of the final model
> 
> cv.sph1<-krige.cv(Residuals~1,data,fit.sph1)
> View(cv.sph1)
> mean(cv.sph1$zscore)  #-0.0009952294
> var(cv.sph1$zscore)   #1.071403
> 
> #Creating a grid file
> 
#==========================================================================
> ==
> 
> grd <- expand.grid(x=seq(from=341960,to=343195,by=10),
>                    y=seq(from=2667253,to=2668873,by=10))
> 
> library(maptools)
> 
> coordinates(grd) <- ~ x+y
> gridded(grd) <- TRUE
> 
> So far,I think, things were going well but I found a error while running
> the following code
> 
> #Regression kriggin
> 
#==========================================================================
> =
> 
> prk.rich<-krige(Richness~Elev, data=data, newdata=grd,model=fit.sph1)
> 
> Error in function (classes, fdef, mtable)  :
>   unable to find an inherited method for function ‘krige’ for signature
> ‘"formula", "missing"’
> 
> Can anyone help me out of this please? I shall be grateful.
> 
> Best regards,
> 
> ---------------
> Md. Abdul Halim
> Assistant Professor
> Department of Forestry and Environmental Science
> Shahjalal University of Science and Technology,Sylhet-3114,
> Bangladesh.
> Cell: +8801714078386.
> alt. e-mail: xou03 at yahoo.com


-----
Ashton Shortridge
Associate Professor			ashton at msu.edu
Dept of Geography			http://www.msu.edu/~ashton
235 Geography Building		ph (517) 432-3561
Michigan State University		fx (517) 432-1671