fit.gstatModel: Methods to fit a regression-kriging model
In GSIF: Global Soil Information Facilities
Description
Usage
Arguments
Details
Note
Author(s)
References
See Also
Examples
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
## S4 method for signature
## '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 signature '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 signature '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 signature '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)
Note
Residuals (response residuals from the model) will be checked for normality and problems reported by default. The warning messages should be taken with care, as when the sample size is small, even big departures from normality will not be reported; when the sample size is large, even the smallest deviation from normality might lead to a warning. Likewise, if the variogram fitting fails, consider fitting a variogram manually or using the fit.vgmModel method.
Author(s)
Tomislav Hengl, Gerard B.M. Heuvelink and Bas Kempen
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.
See Also
gstatModel-class, fit.regModel, test.gstatModel, geosamples-class, stats::glm, gstat::fit.variogram
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# 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)
## End(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)
## End(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))
## End(Not run)
GSIF documentation built on May 2, 2019, 5:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Note Author(s) References See Also Examples
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
## S4 method for signature
## '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 signature '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 signature '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 signature '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 |
formulaString |
object of type |
covariates |
object of type |
method |
character; family of methods considered e.g. |
dimensions |
character; |
fit.family |
character string defyning the GLM family (for more info see |
stepwise |
specifies whether to run step-wise regression on top of GLM to get an optimal subset of predictors |
vgmFun |
variogram function ( |
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 |
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)
Note
Residuals (response residuals from the model) will be checked for normality and problems reported by default. The warning messages should be taken with care, as when the sample size is small, even big departures from normality will not be reported; when the sample size is large, even the smallest deviation from normality might lead to a warning. Likewise, if the variogram fitting fails, consider fitting a variogram manually or using the fit.vgmModel method.
Author(s)
Tomislav Hengl, Gerard B.M. Heuvelink and Bas Kempen
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.
See Also
gstatModel-class, fit.regModel, test.gstatModel, geosamples-class, stats::glm, gstat::fit.variogram
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 | # 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)
## End(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)
## End(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))
## End(Not run)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.