[R-sig-Geo] gstat::krige() - regression kriging vs. kriging with external drift

[Paul Harris]

Hi there
Adding to this…

UK/KED and RK are only equivalent when a global kriging neighbourhood is specified.

If using local kriging neigbourhoods (and assuming quasi-stationarity) then KED/UK can model non-stationary trends (as now the local trend components are filtered out using constraints),
whereas RK doesn’t as it is still using its stationary global trend component in its classic guise (with RK, only the residuals are predicted locally when local neighbourhoods are specified).

However, as noted by Tom, RK is much more flexible in that non-linear trend components can be used.  Median Polish Kriging was probably the first kriging method like this.
KED/UK are always limited to linear trends.

Cheers Harry

-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of rubenfcasal
Sent: 04 December 2015 12:30
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] gstat::krige() - regression kriging vs. kriging with external drift

As I am interested in the topic, I find those comments very useful and I also want to share my thoughts…

 From my point of view UK (Universal Kriging) is a particular case of RK (Regression kriging), UK assumes a linear trend (where spatial coordinates could be used as explanatory variables) . Estimating the trend and computing simple kriging predictions of the residuals is equivalent to UK when a linear trend is assumed (and this trend is estimated efficiently; e.g. Cressie, 1993, section 3.4.5). Note that in most cases ordinary kriging is used instead (e.g if you don't include "
beta = 0" in "gstat::krige(residuals ~ 1, [...])" ).
Regression Kriging is also called Residual Kriging and I would prefer the latter to avoid confusion with other kriging methods that make use of regression techniques.

I don't use the name "Kriging with External Drift" (KED), but I understand that this method considers a linear trend and it is a particular case of UK. Note that the trend could be modelled nonlinearly (or even nonparametrically, see package npsp), so I would preferably say "Kriging with a linear drift".

Note also that kriging methods assume that the variogram is known. In practice, it should be estimated and it is usually done through the residuals, but this is independent of the kriging algorithm…

So I would use (in order of preference): Universal Kriging, Residual Kriging (assuming a linear trend) or Kriging with a linear trend.

Best Regards,

Ruben Fernandez-Casal


El 29/11/2015 a las 20:35, Tomislav Hengl escribió:
>
> I usually use "Kriging with External Drift (KED)" and "Universal
> kriging" as synonyms. I think "Universal kriging" is the oldest term
> (Matheron, 1963), but then 'universal' can seem a bit vague.
>
> Regression-kriging (RK) is when you krige the residuals and add them
> to the trend (hence predict 'm' and 'e' separately). Edzer (gstat)
> implements mainly KED/UK, but then you need the regression residuals
> to estimate the variogram (I guess getting the residuals before you
> fit variogram is unavoidable). Both RK and KED/UK are equivalent i.e.
> give exactly the same predictions under same input settings
> (http://dx.doi.org/10.1016/j.cageo.2007.05.001), so this is where many
> get confused (two paths to the same destination).
>
> The code you have sent is definitively RK.
>
> The real difference between RK and KED is that with RK you can use
> different modelling techniques including random forest and similar. As
> long as the residuals are normal, you can then add interpolated
> residuals back to the trend.
>
> HTH,
>
> T. Hengl
>
> Quoting Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
>
>> On 14/11/15 17:32, Bede-Fazekas Ákos wrote:
>>> Dear List, dear Edzer,
>>> is it correct if I use the term "regression kriging" when I run this
>>> function?:
>>> kriged_value <- gstat::krige(z ~ x + y, [...])@data$var1.pred
>>>
>>> Or should I call it "kriging with external drift" (or "universal
>>> kriging" if x and y are coordinates), and use the term "regression
>>> kriging" only in the case of running this?:
>>> linear_model <- lm(z ~ x + y, [...]) residuals <-
>>> linear_model$residuals kriged_residuals <- gstat::krige(residuals ~
>>> 1, [...])@data$var1.pred kriged_value <- linear_model$fitted.values
>>> + kriged_residuals
>>>
>>
>> For what it's worth, my opinion:
>>
>> I would call it universal kriging, in line with the software you're
>> using (which I wrote):
>>
>>> x = krige(zinc~x+y, meuse, meuse.grid, m)
>> [using universal kriging]
>>
>> Most of the (older) literature associates external drift kriging with
>> a single external drift variable (e.g. Goovaerts), and universal
>> kriging with using coordinates as covariates.
>>
>> The statistician might call everything best linear unbiased
>> prediction under the the general linear model y = X beta + e, with a
>> spatial covariance function describing the covariance of e, and in
>> that perspective both methods are equivalent.
>>
>> I'm not the person to tell what regression kriging is and what it is
>> not, but it might include the above models as well as those where
>> trend fitting and residuals are done in different steps, and
>> potentially under different assumptions for e.
>>
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics (ifgi),  University of Münster,
>> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
>> Journal of Statistical Software:   http://www.jstatsoft.org/
>> Computers & Geosciences:   http://elsevier.com/locate/cageo/
>> Spatial Statistics Society http://www.spatialstatistics.info
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


Rothamsted Research is a company limited by guarantee, registered in England at Harpenden, Hertfordshire, AL5 2JQ under the registration number 2393175 and a not for profit charity number 802038.