GSTAT - measurement vs microscale variation

Edzer,

Thanks again. If both effects (microscale variation and measurement 
error) are collinear in the fit should one not try to include both?

Perhaps it's simply better to jitter by a centimeter. And by the way, 
I'm doing universal kriging with the krige function, is there a way to 
get out the beta coefficients? As in:

krige(formula=no2.ppb~traffic, mydata, newdata=mynewdata,model=vgmRaw)

I'd like to get the WLS or GLS estimate of traffic.

This back and forth is asking a lot of your time, I really appreciate 
your quick replies. Thanks!

Zev

Edzer J. Pebesma wrote:

Zev Ross wrote:

Hi Edzer,

Very useful, thank you. You might be able to tell from my posts that 
I'm running these in parallel in GSTAT and geoR and comparing. It 
seems from your note and example that it might simply be easier to 
add a centimeter (provided a centimeter doesn't matter in the real 
world) to all the coordinates. This way one would not need to keep 
track of the different variances at coincident locations. Do you 
think this would be an acceptable approach?

In terms of your coding, I'm a little uncertain about the meaning of 
the add.to argument and how one might code, for example, half 
micro-scale and half-measurement error. For the example you give, if 
there was a nugget of 0.1 (half micro and half measurement) does my 
coding below look correct? 

Yes

Why is fit.variogram not fitting on the model with error -- it's not 
fitting any of the models I try with measurement error?

Because both effects are collinear in the fit.

I will need to look a bit closer (meaning, in the code), as it 
surprised me that fitting a "Err" effect with no nugget did not work; 
I think it should.

You can find out fit success (to some degree) by looking at the 
"singular" attribute:

myvgmB<-vgm(.5, "Exp",300, add.to=vgm(.05,"Err",0))
fit.variogram(variogram(log(zinc)~1,meuse),model=myvgmB)

 model psill range
1   Err  0.05     0
2   Exp  0.50   300

x = fit.variogram(variogram(log(zinc)~1,meuse),model=myvgmB)
attributes(x)

$names
[1] "model" "psill" "range" "kappa" "ang1"  "ang2"  "ang3"  "anis1" 
"anis2"

$row.names
[1] 1 2

$class
[1] "variogramModel" "data.frame"  $singular
[1] TRUE

$SSErr
[1] 0.0001155682

attr(x, "singular")

[1] TRUE
-- 
Edzer

Zev

myvgmA<-vgm(.5, "Exp",300, nugget=0.1)
myvgmB<-vgm(.5, "Exp",300,nugget=0.05, add.to=vgm(.05,"Err",0))

fit.variogram(variogram(log(zinc)~1,meuse),model=myvgmA)
 model     psill    range
1   Nug 0.0000000   0.0000
2   Exp 0.7186526 449.7581


fit.variogram(variogram(log(zinc)~1,meuse),model=myvgmB)
 model psill range
1   Err  0.05     0
2   Nug  0.05     0
3   Exp  0.50   300





Edzer J. Pebesma wrote:

Zev,

you can use the "Err" variogram model to denote micro variation as 
opposed to nugget. The only effect it has is that for a new 
prediction on an observation location the measurement error-free 
process is predicted, and not the observation process itself. 
Semivariance of an observation with itself remains zero, so it 
doesn't allow for duplicate observations. In terms of predictions, 
it is "as if" you predict right next to a prediction location in 
case the prediction location coincides with an observation location 
(implying that the predicted surface is continuous); in terms of 
prediction variance, it is "as if" you predict for a very small 
block, meaning the nugget is removed from the prediction variance.

Below is an example for the meuse data set.
-- 
Edzer

library(gstat)

Loading required package: sp

data(meuse)
meuse0 = meuse
coordinates(meuse) = ~x+y
# prediction at observation location:
krige(log(zinc)~1,meuse,meuse[1,],vgm(.5, "Exp",300,.5))

[using ordinary kriging]
      coordinates var1.pred var1.var
1 (181072, 333611)  6.929517        0

krige(log(zinc)~1,meuse,meuse[1,],vgm(.5, 

"Exp",300,add.to=vgm(.5,"Err",0)))
[using ordinary kriging]
      coordinates var1.pred  var1.var
1 (181072, 333611)   6.57884 0.1801634

cc = coordinates(meuse)
cc[1,] = cc[1,]+0.01 # 1 cm shift on a 5 km region
coordinates(meuse0)=cc
krige(log(zinc)~1,meuse,meuse[1,],vgm(.5, "Exp",300,.5))

[using ordinary kriging]
      coordinates var1.pred var1.var
1 (181072, 333611)  6.929517        0

krige(log(zinc)~1,meuse,meuse0[1,],vgm(.5, "Exp",300,.5))

[using ordinary kriging]
      coordinates var1.pred var1.var
1 (181072, 333611)  6.578803 0.680188

krige(log(zinc)~1,meuse,meuse0[1,],vgm(.5, 

"Exp",300,add.to=vgm(.5,"Err",0)))
[using ordinary kriging]
      coordinates var1.pred  var1.var
1 (181072, 333611)  6.578803 0.1801880 # same prediction, variance 
0.5 less

krige(log(zinc)~1,meuse,meuse[1,],vgm(.5, 

"Exp",300,.5),block=c(0.01,0.01))
[using ordinary kriging]
      coordinates var1.pred  var1.var
1 (181072, 333611)  6.578836 0.1801594



Zev Ross wrote:

Hi All,

I folded this question into a previous post, but I think it may 
have gotten
missed. Just wondering if someone could tell me how, in GSTAT, one 
would specify
the nugget as measurement error vs microscale variation in 
kriging. I have
multiple measurements at the same location and I'd like to use 
these to
determine the measurement error. I've figured out how to do this 
in geoR, but as
most of my scripts are written in R using GSTAT, I'd rather use that.

Thank you! Zev