[R-sig-Geo] indicator kriging with R

[Edzer Pebesma]

On 23/10/15 01:45, Tim Peterson wrote:
> Hi Dafne and Edzer,
> 
> I too have been looking into gstat indicator kriging and simulation. My 
> interest is because indicator simulations can (at least in GSLib) 
> account for systems where extreme values are more/less spatially 
> correlated than median values; that is the upper and lower threshold 
> variograms can have a very different range to the mean values.
> 
> However, following 
> http://rstudio-pubs-static.s3.amazonaws.com/10216_57cb0ffbb2e94586b29dfb1849dd49f0.html, 
> I've just realised that gstat indicator kriging requires all indicator 
> variograms to have an intrinsic correlation or be a linear model of 
> coregionalisation, and hence requiring each variogram to have at least 
> an equal range (FYI, to test this change line 'meuse.fit = fit.lmc(x, 
> meuse.i)' to 'meuse.fit = fit.lmc(x, meuse.i,fit.ranges=TRUE)'). This is 
> most likely because, unlike GSLib, gstat uses cokriging but this 
> limitation defeats one of the advantages of indicator simulation over 
> Gaussian simulation.

please look into the docs of ?fit.lmc, which lets you fit non-LMC
models, with each variogram its own (fitted) range. You can use these
models for direct kriging, and also for cokriging, but in the latter
case there is of course no guarantee that the covariance matrix is PD -
the software will issue warnings, you're at your own risk, and this is
quite likely hard to defend for reviewers who know some basics of random
processes.

> 
> Anyway, can anyone propose a work around to this limitation? The only 
> solution I can see is to krige each threshold independently but this 
> does eliminate the option of accounting for cross-correlations between 
> the thresholds or more advanced indicator feature; such as including 
> spatially varying constraints eg 
> http://www.mssanz.org.au/modsim2011/I9/peterson.pdf

Another, more consistent approach to modeling spatial dependencies where
correlation strength varies across the value range uses copula's, see
e.g. Ben Graeler's PhD work, http://ifgi.uni-muenster.de/~b_grae02/ and
http://spcopula.r-forge.r-project.org/

> 
> Cheers
> 
> Tim
> 
> 
> On 22/10/15 20:12, Dafne Berg wrote:
>> Dear Edzer,
>> yes it is exactly what I meant. Thanks for the confirmation that it is not
>> yet implemented in gstat. I will have a look to gslib and the book.
>>
>> Thanks a lot
>>
>> Dafne
>>
>> On 22 October 2015 at 09:38, Edzer Pebesma <edzer.pebesma at uni-muenster.de>
>> wrote:
>>
>>> Yes, this is what indicator kriging does.
>>>
>>> What you want is probably going back from the [0..1] estimated
>>> probabilities to the original Z values. Gstat doesn't help here, but it
>>> shouldn't be too hard to write a function for this. You need to make
>>> assumptions about the distribution between different cutoff values, in
>>> particular about both extreme (open) classes. I think the GSLIB and
>>> Goovaerts books would be the first ones to look into.
>>>
>>> On 22/10/15 10:19, Dafne Berg wrote:
>>>> Dear list,
>>>> I am trying to do a indicator kriging of a continuous variable. I am
>>>> especially interested in performing the "order relation correction". In
>>> the
>>>> documentation of "predict" in gstat it is mentioned as debug level 64
>>>> (order relation violations (indicator kriging values before and after
>>> order
>>>> relation correction).
>>>>
>>>> However, when running predict I can only get the values in the range 0 to
>>>> 1, as expected for indicator kriging.
>>>>
>>>> I am probably missing something very basic here, but I will be grateful
>>> for
>>>> any pointers in the right direction
>>>>
>>>> Thank you very much in advance
>>>>
>>>> Best wishes
>>>>
>>>> Dafne
>>>>
>>>>        [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>> --
>>> 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
>>>
>> 	[[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151023/4a6f9415/attachment.bin>