1

I wish to estimate a variogram for spatially distributed price data in Tanzania. I'm new to both spatial statistics and gstat, and have a couple questions.

(1) I can SEE spatial patterns in price, and moran's I is about 0.2, so I know there is spatial correlation. But when I plot my variogram using the code below, semivariance is flat across distance. (Screenshot in dropbox link below, as well as the tzprice1 shapefile for replication.) Why would this be? Perhaps I need to specify smaller intervals, and a smaller range? And if so, how?

tzprice1_v <- tzprice1[!is.na(tzprice1@data$prodp), ]
pr.v <- variogram(prodp ~1, tzprice1_v)
plot(pr.v)

(2) I am unable to estimate the variogram using fit.variogram, and I'm guessing that this is because the underlying data don't show the right shape. Is this correct? But just in case, I'm using the code below... is it correct to simply guess the sill, nugget, and range starting values in this way, from the plot that I did above? 

fit.variogram(pr.v, vgm(psill=.125,"Exp",range=200,nugget=.1))

Shapefile and screenshot here: https://www.dropbox.com/sh/oad1xfcpcugkulh/AACec_FuUZqCHR9-Jicohgm9a?dl=0

| improve this question | |
5

If you make your variogram go a bit further you can start to see some structure:

> pr.v <- variogram(prodp ~1, tzprice1_v, cutoff=1500, width=1500/8)
> plot(pr.v)

enter image description here

You can then fit a variogram if you use a different fit weighting method:

> fv = fit.variogram(pr.v, vgm(psill=.05,"Exp",range=80,nugget=.15),fit.method=1)
> plot(pr.v, fv)

enter image description here

Or with fit.method=6:

enter image description here

which seems to give the last point enough weight to bring the curve down a bit.

The next thing to do is to compute your predictions and see if either of these variograms make any significant difference to the predictions and SEs.

What is going on is that you have spatial correlation over a large scale which seems to be due to trends in the data - you might do better to fit a linear or polynomial trend surface when kriging to soak up most of the large-scale correlation.

| improve this answer | |
  • In the fit.variogram() function vgm(psill=.05,"Exp",range=80,nugget=.15) looks a little bit strange. (I mean, too small psill and too small range.) A typo, maybe? – Kazuhito Mar 20 '18 at 1:10
  • No typo, that's the initial value, the fitted value looks to be very different - I should maybe have printed it out. – Spacedman Mar 20 '18 at 8:01
  • I agree we do not have to worry too much about the initial value in vgm(). Then vgm("Exp") could be enough?... if we are using the latest gstat. – Kazuhito Mar 20 '18 at 8:40
  • This is incredibly, incredibly helpful. Thank you SO much! A few quick follow-ups. (1) This distance is in km, correct? Since unprojected? (2) If so, the semivariance data -- never mind the fitted function -- implies that the bulk of spatial correlation happens w/in 900-1000km, right? That's where the semiv. plateaus. (3) But the fitted model gives a range of ~519. Why? (4) The fitted model seems poor, to me; a polynomial would allow the slope to descend. Why are the only options linear and exponential? – Leah Bevis Mar 23 '18 at 19:50
1

I would go for the automap package:

library(automap)
fit.vgm = autofitVariogram(prodp ~1, tzprice1_v, model = "Exp")
plot(pr.v, fit.vgm$var_model)

enter image description here

On the other hand though, I'm not a statistician - I work with variograms more empirically which I guess I shouldn't recommend. I am not able to answer your question about possible ways of estimating variogram parameters, nor am I able to tell which variogram model would me the most appropriate (it seems like some sort of Wave pattern in your data but this theoretical variogram didn't fit to your empirical variogram).

| improve this answer | |
  • In a Kriging model, this semivariogram would be estimating nothing but error (pure nugget effect). Please see the post from Spacedman on some semivariograms that, given the data, would be providing more correct estimates. – Jeffrey Evans Mar 19 '18 at 15:34
  • @JeffreyEvans Thanks - yes, those variograms look much better :-) As I said, I'm using it more mechanically as I employ kriging in routine operations where I can't manually check the data. This was just one of the options that R offers. – Janina Mar 20 '18 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.