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

[Zev Ross]

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

-- 
Zev Ross
ZevRoss Spatial Analysis
303 Fairmount Ave
Ithaca, NY 14850
607-277-0004 (phone)
866-877-3690 (fax, toll-free)
zev at zevross.com