gstat now uses Lapack / failing cokriging

Hi Andrew and Edzer,

I checked the gstat manual: on page 70 it writes that it uses the 'classic' normal cross-variogram when the two variables have the same number of observations and identical coordinates and order, and that in all other cases it uses the pseudo cross-variogram. In the example that Sytze and I used it must then have used the normal cross-variogram.

Personally I am not in favour of using the pseudo cross-variogram because it messes up the measurement units (it is defined as Gij(sk-sl) = 0.5*E[(Zi(sk) - Zj(sl))^2], where for example Zi might be rainfall in mm and Zj temperature in Kelvin), but I also see the problems with the normal cross-variogram in that it cannot reproduce non-symmetric cross-covariance functions and can only be estimated from paired observations at locations where both variables were measured. Note that the suboptimal performance when using the normal cross-variogram as reported in Ver Hoef and Cressie only occurs when the cross-covariance function is non-symmetric (i.e. when in their example parameter DELTA is not equal to zero).

When gstat estimates the pseudo cross-variogram then I think that it should derive from it the cross-covariances as Cij(sk-sl) = 0.5*Cii(0) + 0.5*Cjj(0) - Gij(sk-sl) (Ver Hoef and Cressie page 234), i.e. it should use the average of the sill values of the two variograms and not max(0, sum of the positive sill values).

Gerard

-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of r-sig-geo-request at r-project.org
Sent: 14 December 2015 12:00
To: r-sig-geo at r-project.org
Subject: R-sig-Geo Digest, Vol 148, Issue 12

Send R-sig-Geo mailing list submissions to
	r-sig-geo at r-project.org

To subscribe or unsubscribe via the World Wide Web, visit
	https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
	r-sig-geo-request at r-project.org

You can reach the person managing the list at
	r-sig-geo-owner at r-project.org

When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-Geo digest..."


Today's Topics:

   1. "Open Source software and web-based apps for soil data",
      session at EGU 2016 (Tomislav Hengl)
   2. Re: gstat now uses Lapack / failing cokriging
      (Andrew Zammit Mangion)
   3. Re: FW: gstat now uses Lapack / failing cokriging (Edzer Pebesma)


------------------------------

Message: 2
Date: Mon, 14 Dec 2015 00:18:53 +0000
From: Andrew Zammit Mangion <azm at uow.edu.au>
To: "'r-sig-geo at r-project.org'" <r-sig-geo at r-project.org>
Subject: Re: [R-sig-Geo] gstat now uses Lapack / failing cokriging
Message-ID: <1450052325363.12086 at uow.edu.au>
Content-Type: text/plain; charset="iso-8859-1"

Hi Gerard,
  Just a small comment on the use of pseudo- vs. normal- cross-semivariograms. Normal cross-semivariograms may yield nonoptimal weights for cokriging, as shown in the example of Ver Hoef and Cressie (1993), and I believe it's a nomenclature problem that the 'pseudo' one is that which should be used. A quick look at the gstat documentation I see that the pseudo cross-semivariogram is what is used but maybe Edzer can confirm.

Best,
Andrew

-----------------------------------------------------------------------------
Andrew Zammit Mangion, 
School of Mathematics and Applied Statistics, 
University of Wollongong NSW 2522
------------------------------------------------------------------------------

Hi Edzer,

I checked the Ver Hoef and Cressie (1993) paper but could not find a reference to deriving the cross-covariance from the cross-variogram as "C(h)=C(0)-gamma(h) with C(0)=max(0, sum of the positive sill values)". I did notice that their paper is mostly about the pseudo cross-variogram (defined in Eq. 9), while we and presumably gstat as well work with the 'normal' cross-variogram (Eq. 8).

I think the proper relation is Cij(h)=Cij(0)-Vij(h), where Cij(h)=E[(Zi(s)-E[Zi(s)])*(Zj(s+h)-E[Zj(s+h)])] is the cross-covariance and Vij(h)=0.5*E[(Zi(s)-Zi(s+h))*(Zj(s)-Zj(s+h))] is the normal cross-variogram. This relation is valid if Z is second-order stationary (conditions 1 and 2 on page 220 of VH&C) and satisfies the symmetric cross-covariance condition (i.e. Cij(h)=Cij(-h), condition 3 on the same page). I would recommend that you implement this relation instead of "C(h)=C(0)-gamma(h) with C(0)=max(0, sum of the positive sill values)".

What I still do not understand is that the 'old' version of gstat did not complain about our Allier example, because using the alternative derivation of the cross-covariance from the cross-variogram produces a kriging system that is not positive-definite.

Gerard