[R-sig-Geo] fit.lmc in gstat - fit same range

Edzer Pebesma edzer.pebesma at uni-muenster.de
Sun Nov 23 12:34:09 CET 2008


Hi Ole,

You're completely right -- it is hard to explain why the model you get 
when fitting different ranges is useful. If you had tried to use it in a 
kriging context, you would have seen the following:

 > predict(meuse.fit3, meuse[1,])
Warning: No Intrinsic Correlation or Linear Model of Coregionalization found
Reason: ranges differ
Warning: [add `set nocheck = 1;' to the command file to ignore the 
following error]

Error in predict.gstat(meuse.fit3, meuse[1, ]) :
  gstat: value not allowed for: variograms do not satisfy a legal model

Only few on this list will know that you can override this error by:
 > meuse.fit3$set=list(nocheck=1)
 > predict(meuse.fit3, meuse[1,])
Warning: No Intrinsic Correlation or Linear Model of Coregionalization found
Reason: ranges differ
Now checking for Cauchy-Schwartz inequalities:
variogram(var0,var1) passed Cauchy-Schwartz
variogram(var0,var2) passed Cauchy-Schwartz
variogram(var1,var2) passed Cauchy-Schwartz
variogram(var0,var3) passed Cauchy-Schwartz
variogram(var1,var3) passed Cauchy-Schwartz
variogram(var2,var3) passed Cauchy-Schwartz
[using ordinary cokriging]
       x      y  zn.pred       zn.var  cu.pred       cu.var  cd.pred
1 181072 333611 6.929517 1.004906e-26 4.442651 1.665335e-16 2.459589
        cd.var  pb.pred       pb.var    cov.zn.cu    cov.zn.cd    cov.cu.cd
1 8.881784e-16 5.700444 1.110223e-16 1.110223e-16 4.440892e-16 3.330669e-16
     cov.zn.pb    cov.cu.pb    cov.cd.pb
1 1.110223e-16 5.551115e-17 1.110223e-15

which may be useful to some. IIRC, I left this option in because some 
people at some stage advocated to use cokriging of indicators for 
several thresholds of a single continuous variable, where the LMC makes 
less sense because extremes show smaller correlation distance--all this 
is definitely not in the Webster & Oliver book.

The reason why fit.lmc was not cleanly implemented the way you would 
like is that it uses consecutive calls to fits for individual 
variograms, and if necessary corrects for positive definiteness of 
sills/nuggets. This setup doesn't make fitting of a common sill easy. 
Nobody has volunteered to do it better so far.

A larger worry I have now is that the gstat code allows you to predict 
values on a sphere with all the regular variogram models that are not 
positive on the sphere--without this kind of error/warning.
--
Edzer

Ole Fredslund Christensen wrote:
> Dear experts
>
> I'm in the middle of giving a course on geostatistics and GIS (I am doing the geostatistics part, 
> someone else the GIS part). We are using Webster and Oliver as text book.
>
> I am considering to make a small exercise related to Multivariate geostatistics and cokriging.
> In Webster and Oliver they present the linear model of coregionalisation 
> such that 
> 1) it is a sum of components
> 2) for component the covariances and crosscovarainces share same correlation function,
> but have different scalars associated to them (these scalar are chosen consistently to make 
> up a covariance matrix) I can understand why this construction is a valid one 
> (and I hope to able to explain it to the students)
>
>
> In gstat I can do the following :
>
> data(meuse)
>
> pdf(file="meuse.pdf")
> meuse.g <- gstat(id = "zn", formula = log(zinc) ~ 1, loc = ~x + y, data = meuse)
> meuse.g <- gstat(meuse.g, "cu", log(copper) ~ 1, ~x + y, meuse)
> meuse.g <- gstat(meuse.g, "cd", log(cadmium) ~ 1, ~x + y, meuse)
> meuse.g <- gstat(meuse.g, "pb", log(lead) ~ 1, ~x + y, meuse)
>
> data(meuse)
>
> pdf(file="meuse.pdf")
> meuse.g <- gstat(id = "zn", formula = log(zinc) ~ 1, loc = ~x + y, data = meuse)
> meuse.g <- gstat(meuse.g, "cu", log(copper) ~ 1, ~x + y, meuse)
> meuse.g <- gstat(meuse.g, "cd", log(cadmium) ~ 1, ~x + y, meuse)
> meuse.g <- gstat(meuse.g, "pb", log(lead) ~ 1, ~x + y, meuse)
>
>
> meuse.g2 <- gstat(meuse.g, model = vgm(1, "Sph", 900,1), fill.all = T)
>   
you forgot:
x = variogram(meuse.g)
> meuse.fit2 <- fit.lmc(x, meuse.g2)
> meuse.fit2
>
> variograms:
>          model      psill range
> zn[1]      Nug 0.05838220     0
> zn[2]      Sph 0.57717809   900
> cu[1]      Nug 0.06127586     0
> cu[2]      Sph 0.24767761   900
> cd[1]      Nug 0.48904863     0
> cd[2]      Sph 1.22387180   900
> pb[1]      Nug 0.05924609     0
> pb[2]      Sph 0.47945611   900
> zn.cu[1]   Nug 0.04630670     0
> zn.cu[2]   Sph 0.35863043   900
> zn.cd[1]   Nug 0.12881136     0
> zn.cd[2]   Sph 0.78603818   900
> cu.cd[1]   Nug 0.12206442     0
> cu.cd[2]   Sph 0.49782101   900
> zn.pb[1]   Nug 0.05197207     0
> zn.pb[2]   Sph 0.51367619   900
> cu.pb[1]   Nug 0.03932104     0
> cu.pb[2]   Sph 0.30882315   900
> cd.pb[1]   Nug 0.10175399     0
> cd.pb[2]   Sph 0.68735760   900
> ~x + y
>
>
> ## same range, but it not fit.
>
> meuse.fit3 <- fit.lmc(x, meuse.g2, fit.ranges =TRUE)
> mesure.fit3
>
> variograms:
>          model      psill     range
> zn[1]      Nug 0.06390664    0.0000
> zn[2]      Sph 0.59513904  965.6013
> cu[1]      Nug 0.05466579    0.0000
> cu[2]      Sph 0.23630293  773.4431
> cd[1]      Nug 0.51350775    0.0000
> cd[2]      Sph 1.32830734 1063.0889
> pb[1]      Nug 0.07325686    0.0000
> pb[2]      Sph 0.56031923 1193.4988
> zn.cu[1]   Nug 0.04619022    0.0000
> zn.cu[2]   Sph 0.35832564  898.0459
> zn.cd[1]   Nug 0.14318935    0.0000
> zn.cd[2]   Sph 0.84453467 1044.6003
> cu.cd[1]   Nug 0.12402376    0.0000
> cu.cd[2]   Sph 0.50344066  924.9874
> zn.pb[1]   Nug 0.06218012    0.0000
> zn.pb[2]   Sph 0.55728401 1062.1928
> cu.pb[1]   Nug 0.04216347    0.0000
> cu.pb[2]   Sph 0.31801824  962.8714
> cd.pb[1]   Nug 0.11849857    0.0000
> cd.pb[2]   Sph 0.76933704 1117.9877
>
>
> fit different ranges !!!
> (and I have difficulties to understand why this construction neccessarily is valid, 
> and clearly would have problems explaining it to students).
>
> What I really would like to be able to do is to fit the range also, but restrict it to be the same for
> all varoigrams and crossvariograms
> (requring a simultaneous fitting of variograms and cross-variograms, I guess).
>
> I can see a solution,  where I chose 5 different fixed ranges and make a visual inspection of this
> for each of them, but I worry about draging students through an exercise where they inspect a number of 
> different variograms and for each of these have to consider five fitted variograms
>
> Anyone aware of a more elegant solution to this problem ?
>
>
> Ole Christensen
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de




More information about the R-sig-Geo mailing list