[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