[R-sig-Geo] GSTAT - measurement vs microscale variation

[Edzer J. Pebesma]

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
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>   
>>
>>
>>
>