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

[Ole Fredslund Christensen]

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