view raw Rmd

Fitting variogram functions with R package gstat has become more flexible, and hopefully more user friendly. Up to now, after loading data

library(sp)
demo(meuse, ask = FALSE, echo = FALSE) # load meuse data set

users were required to use a sequence like

library(gstat)
v = variogram(log(zinc)~1, meuse)
v.fit = fit.variogram(v, vgm(1, "Sph", 900, 1))
v.fit

##   model      psill    range
## 1   Nug 0.05066243   0.0000
## 2   Sph 0.59060780 897.0209

where fit.variogram fits variogram parameters of a spherical model (Sph) to the sample variogram v. The values 1, 900 and 1 were needed as initial values in the weighted non-linear fit (where only the range parameter is non-linear).

This has changed in gstat version 1.2: now, vgm can take only a variogram model, as in

fit.variogram(v, vgm("Sph"))

##   model      psill    range
## 1   Nug 0.05065971   0.0000
## 2   Sph 0.59060511 897.0011

or even a set of models, in which case the best fitting is returned, as in

fit.variogram(v, vgm(c("Exp", "Sph")))

##   model      psill    range
## 1   Nug 0.05065971   0.0000
## 2   Sph 0.59060511 897.0011

fit.variogram(v, vgm(c("Exp", "Mat", "Sph")))

##   model      psill    range
## 1   Nug 0.05065971   0.0000
## 2   Sph 0.59060511 897.0011

where we still see that the sperical model is chosen. If we choose a different sample variogram, where Matern is chosen, as in:

v0 = variogram(zinc~1, meuse)
fit.variogram(v0, vgm(c("Exp", "Mat", "Sph")))

##   model    psill    range kappa
## 1   Nug   9486.4   0.0000   0.0
## 2   Mat 163285.3 381.7076   0.5

we see that the kappa value is 0.5, which is a default value that was not fit. We can fit kappa by specifying fit.kappa = TRUE, as in

options(warn = -1) # don't print warnings
fit.variogram(v0, vgm(c("Exp", "Mat", "Sph")), fit.kappa = TRUE)

##   model    psill    range kappa
## 1   Nug      0.0   0.0000   0.0
## 2   Mat 176455.5 456.4375   0.4

where the best fitting kappa from the range 0.3, 0.4, 0.5,…,5 is chosen. I suppressed warnings here, as around 20 warnings were printed in cases with crazy initial values. This is usual for Matern models: larger kappa values have effective ranges (the distance value at which the model reaches, say, 95% of its sill) much larger than the range parameter, as illustrated by

plot(variogramLine(vgm(1, "Mat", 1, kappa = 4), 10), type = 'l')

where at distance 1, 0.05 of the sill is reached (and the model, up till there, is nearly linear or parabolic, leading to a singularity during fit). A different parameterisation of the Matern model, given in Michael Stein’s book, is the following

plot(variogramLine(vgm(1, "Ste", 1, kappa = 4), 10), type = 'l')

This one has the same smoothness, but reaches the sill much closer to the range value. As a consequence it fits easier, that is, without warnings:

options(warn = 0) # set back to normal
fit.variogram(v0, vgm(c("Exp", "Ste", "Sph")), fit.kappa = TRUE)

##   model    psill    range kappa
## 1   Nug      0.0   0.0000   0.0
## 2   Ste 176455.4 577.3522   0.4

For those you need a more precise estimate of the optimal kappa value, you can iterate over steps of e.g. 0.01 by

fit.variogram(v0, vgm(c("Exp", "Ste", "Sph")), fit.kappa = seq(.3,5,.01))

##   model       psill    range kappa
## 1   Nug    282.7529   0.0000  0.00
## 2   Ste 175030.7434 563.3126  0.41

How it works

Default initial parameter values are chosen from the sample variogram, where:

  • the range parameter is taken as 1/3 of the maximum sample variogram distance,
  • the nugget parameter is taken as the mean of the first three sample variogram values, and
  • the partial sill is taken as the mean of the last five sample variogram values.
vgm("Sph")

##   model psill range
## 1   Nug    NA     0
## 2   Sph    NA    NA

contains NA values for the numeric parameters, and under the hood (undocumented)

gstat:::vgm_fill_na(vgm("Sph"), v)

##   model     psill    range
## 1   Nug 0.2141508   0.0000
## 2   Sph 0.6200691 514.4008

fills the NA values with the initial values.

Providing more than one model to vgm returns a list,

vgm(c("Sph", "Exp"))

## [[1]]
##   model psill range
## 1   Nug    NA     0
## 2   Sph    NA    NA
## 
## [[2]]
##   model psill range
## 1   Nug    NA     0
## 2   Exp    NA    NA
## 
## attr(,"class")
## [1] "variogramModelList" "list"

which fit.variogram iterates over, returning the best fitting model.

Comparison to automap

Function automap::autofitVariogram does a similar job, but includes the computation of the sample variogram from data (which can be controlled by passing parameters to ...). It takes slightly different defaults for fitting, definitely different defaults when computing the sample variogram, and has options for combining distance bins.