fit.gstatModel-methods: Methods to fit a regression-kriging model

Description

Tries to automatically fit a 2D or 3D regression-kriging model for a given set of points (object of type "SpatialPointsDataFrame" or "geosamples") and covariates (object of type "SpatialPixelsDataFrame"). It first fits a regression model (e.g. Generalized Linear Model, regression tree, random forest model or similar) following the formulaString, then fits variogram for residuals usign the fit.variogram method from the gstat package. Creates an output object of class gstatModel-class.

Usage

# S4 method for SpatialPointsDataFrame,formula,SpatialPixelsDataFrame
fit.gstatModel(observations, formulaString, covariates, 
     method = list("GLM", "rpart", "randomForest", "quantregForest",
      "xgboost", "ranger"), 
     dimensions = list("2D", "3D", "2D+T", "3D+T"),
     fit.family = gaussian(), stepwise = TRUE, vgmFun = "Exp", 
     subsample = 5000, subsample.reg = 10000, …)
# S4 method for geosamples,formula,SpatialPixelsDataFrame
fit.gstatModel(observations, formulaString, covariates, 
     method = list("GLM", "rpart", "randomForest", "quantregForest", 
     "xgboost", "ranger"), 
     dimensions = list("2D", "3D", "2D+T", "3D+T"),
     fit.family = gaussian(), stepwise = TRUE, 
     vgmFun = "Exp", subsample = 5000, subsample.reg = 10000, …)
# S4 method for geosamples,formula,list
fit.gstatModel(observations, formulaString, covariates, 
     method = list("GLM", "rpart", "randomForest", "quantregForest", 
     "xgboost", "ranger"), 
     dimensions = list("2D", "3D", "2D+T", "3D+T"),
     fit.family = gaussian(), stepwise = TRUE, 
     vgmFun = "Exp", subsample = 5000, subsample.reg = 10000, …)
# S4 method for geosamples,list,list
fit.gstatModel(observations, formulaString, covariates, 
     method = list("GLM", "rpart", "randomForest", "quantregForest",
      "xgboost", "ranger"), 
     dimensions = list("2D", "3D", "2D+T", "3D+T"),
     fit.family = gaussian(), stepwise = TRUE, 
     vgmFun = "Exp", subsample = 5000, subsample.reg = 10000, …)

Arguments

observations

object of type "SpatialPointsDataFrame" or "geosamples-class"

formulaString

object of type "formula" or a list of formulas

covariates

object of type "SpatialPixelsDataFrame", or list of grids

method

character; family of methods considered e.g. "GLM"

dimensions

character; "3D", "2D", "2D+T", "3D+T" models

fit.family

character string defyning the GLM family (for more info see stats::glm)

stepwise

specifies whether to run step-wise regression on top of GLM to get an optimal subset of predictors

vgmFun

variogram function ("Exp" by default)

subsample

integer; maximum number of observations to be taken for variogram model fitting (to speed up variogram fitting)

subsample.reg

integer; maximum number of observations to be taken for regression model fitting (currently only used for randomForest)

…

other optional arguments that can be passed to glm and/or fit.variogram

Details

The GLM method by default assumes that the target variable follows a normal distribution fit.family = gaussian(). Other possible families are:

normal distribution: fit.family = gaussian() (default setting)
log-normal distribution: fit.family = gaussian(log)
binomial variable: fit.family = binomial(logit)
variable following a poisson distribution: fit.family = poisson(log)

References

Meinshausen, N. (2006). Quantile regression forests. The Journal of Machine Learning Research, 7, 983-999.
chapter 8 ``Interpolation and Geostatistics'' in Bivand, R., Pebesma, E., Rubio, V., (2008) Applied Spatial Data Analysis with R. Use R Series, Springer, Heidelberg, pp. 378.
Hengl, T. (2009) A Practical Guide to Geostatistical Mapping, 2nd Edt. University of Amsterdam, www.lulu.com, 291 p.

Examples

Run this code

# NOT RUN {
# 2D model:
library(sp)
library(boot)
library(aqp)
library(plyr)
library(rpart)
library(splines)
library(gstat)
library(randomForest)
library(quantregForest)
library(plotKML)

## load the Meuse data set:
demo(meuse, echo=FALSE)

## simple model:
omm <- fit.gstatModel(meuse, om~dist+ffreq, meuse.grid,
   family = gaussian(log))
om.rk <- predict(omm, meuse.grid)
plot(om.rk)
## it was succesful!

## fit a GLM with a gaussian log-link:
omm <- fit.gstatModel(meuse, om~dist+ffreq, meuse.grid, 
   fit.family = gaussian(log))
summary(omm@regModel)
om.rk <- predict(omm, meuse.grid)
plot(om.rk)

## fit a regression-tree:
omm <- fit.gstatModel(meuse, log1p(om)~dist+ffreq, meuse.grid, 
   method="rpart")
summary(omm@regModel)
## plot a regression-tree:
plot(omm@regModel, uniform=TRUE)
text(omm@regModel, use.n=TRUE, all=TRUE, cex=.8)
omm@vgmModel    

## fit a randomForest model:
omm <- fit.gstatModel(meuse, om~dist+ffreq, meuse.grid, 
   method="randomForest")
## plot to see how good is the fit:
plot(omm)
## plot the estimated error for number of bootstrapped trees:
plot(omm@regModel)
omm@vgmModel
om.rk <- predict(omm, meuse.grid)
plot(om.rk)
## Compare with "quantregForest" package:
omm <- fit.gstatModel(meuse, om~dist+ffreq, meuse.grid, 
   method="quantregForest")
# }
# NOT RUN {
om.rk <- predict(omm, meuse.grid, nfold=0)
plot(om.rk)
## plot the results in Google Earth:
plotKML(om.rk)
# }
# NOT RUN {
## binary variable (0/1):
meuse$soil.1 <- as.numeric(I(meuse$soil==1))
som <- fit.gstatModel(meuse, soil.1~dist+ffreq, meuse.grid, 
   fit.family = binomial(logit))
summary(som@regModel)
som.rk <- predict(som, meuse.grid)
plot(som.rk)
# }
# NOT RUN {
# plot the results in Google Earth:
plotKML(som.rk)
# }
# NOT RUN {
## 3D model:
library(plotKML)
data(eberg)
## list columns of interest:
s.lst <- c("ID", "soiltype", "TAXGRSC", "X", "Y")
h.lst <- c("UHDICM","LHDICM","SNDMHT","SLTMHT","CLYMHT")
sel <- runif(nrow(eberg))<.05
## get sites table:
sites <- eberg[sel,s.lst]
## get horizons table:
horizons <- getHorizons(eberg[sel,], idcol="ID", sel=h.lst)
## create object of type "SoilProfileCollection"
eberg.spc <- join(horizons, sites, type='inner')
depths(eberg.spc) <- ID ~ UHDICM + LHDICM
site(eberg.spc) <- as.formula(paste("~", paste(s.lst[-1], collapse="+"), sep=""))
coordinates(eberg.spc) <- ~X+Y
proj4string(eberg.spc) <- CRS("+init=epsg:31467")
## convert to geosamples:
eberg.geo <- as.geosamples(eberg.spc)
## covariates:
data(eberg_grid)
gridded(eberg_grid) <- ~x+y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")
glm.formulaString = as.formula(paste("SNDMHT ~ ", 
  paste(names(eberg_grid), collapse="+"), "+ ns(altitude, df=4)"))
SNDMHT.m <- fit.gstatModel(observations=eberg.geo, glm.formulaString, 
  covariates=eberg_grid)
plot(SNDMHT.m)
## problems with the variogram?
# }
# NOT RUN {
## remove classes from the PRMGEO6 that are not represented in the model:
sel = !(levels(eberg_grid$PRMGEO6) %in% levels(SNDMHT.m@regModel$model$PRMGEO6))
fix.c = levels(eberg_grid$PRMGEO6)[sel]
summary(eberg_grid$PRMGEO6)
for(j in fix.c){
  eberg_grid$PRMGEO6[eberg_grid$PRMGEO6 == j] <- levels(eberg_grid$PRMGEO6)[7]
}
## prepare new locations:
new3D <- sp3D(eberg_grid)
## regression only:
SNDMHT.rk.sd1 <- predict(SNDMHT.m, new3D[[1]], vgmmodel=NULL)
## regression-kriging:
SNDMHT.rk.sd1 <- predict(SNDMHT.m, new3D[[1]])
## plot the results in Google Earth:
plotKML(SNDMHT.rk.sd1, z.lim=c(5,85))
# }

Run the code above in your browser using DataLab