An exponential spatial function was fitted to experimental variogram of normalized yield variable.The result obtained look like:
## Exponential variogram model
> exp.var exp.mod print(exp.mod)
model psill range
1 Nug 0.0000000 0.00000
2 Exp 0.8066099 11.36279
The function does not exhibit nugget effect. That is, the nugget = 0. Does it mean that the exponential structure accounted for 100 of spatial variation? The attached file shows the semivariogram plot fitted with exponential function. Is there a way of showing the horizontal line in gstat to indicate theoretical sill? It seems the data point rises without level off. Does it mean there is a trend in the data? Please have a look at the file and comment.
Thanks
It is difficult to tell if the variogram fits the data just looking at the empirical semi variogram plot. It may be that there are few pairs of data for m>100 but that is not clear in the plot. You can fix the sill value and rerun the data to see if the plots looks better, or try another model.
You can change the fit method, for example, fit.variogram(fit.method=7,...) is the default method. The other methods are given here:
See http://www.gstat.org/gstat.pdf table 2 for more information.
fit by weight
0 - - (no t)
1 gstat Nj
2 gstat Nj/gamma(hj)^2
3 gnuplot Nj
4 gnuplot Nj/gamma(hj)^22
5 gstat REML
6 gstat no weights (OLS)
7 gstat Nj/h^2
If you want to fit this plot use fit.method=6,
There are other packages that I would try as well here:
Such as geoR:
http://www.leg.ufpr.br/geoR/geoRdoc/geoRintro.html
The plots generate trellis obejcts. You may be able to use the latticeExtra package to add lines to this plot, otherwise you have to modify the panel functions defined by
gstat:::plot.gstatVariogram
Looks like a poor fit to me. I think two structures would fit these estimates almost exactly. (You may also do better with a single structure but with a much longer range of influence) Regarding a trend, the lack of perfectly horizontal sill does indicate some larger scale gradient. I would do a regression on x and y, or other covariates you may have a available. That said, the shallow slope of the increasing sill just suggests to me that the variance explained by any trend you fit will be minimal and probably unimportant if you are going to use this variogarm for Kriging.
Looks to me like your fitting routine is finding a local minimum, so I would experiment with a range of starting guesses for the fitting routine. I would also experiment with selecting some parameters by hand and plotting to get a feel for how close the Exp model can come to fitting both the long range and short range behavior of this curve.
Fist of all the exponential model only has an asymptotic sill. Secondly, when you compute the experimental variogram ("semivariogram" is long outdated terminology) it is important to only compute it for a total distance of about one half the diagonal distance across the region determined by the data locations. Beyond that distance the number of pairs begins to decline rapidly and the plotted points are not reliable estimates of the variogram values.
Before computing the experimental variogram, try fitting a trend surface to the data both for degree one polynomials and for second degree polynomials. If the coefficients appear to be non-zero then try using the residuals to compute the experimental variogram and compare with the plot you already have. If you have not already done so, make a coded plot of the data locations, each data location coded by the data value at that point. This is a visual way to check for a non-constant drift. The experimental variogram only estimates the values of the variogram IF the expected value of the random function is a constant. The trend is not theoretically the same as the drift but is commonly used to estimate it. If the experimental variogram appears to grow at a quadratic or higher rate that is an indication of a non-constant drift..
Note that the sample variance is not a good estimate of the sill, it is based on assumption of independence.
I would also suggest using the model you got from gstat and look at the cross-validation results,e.g, is the arithmetic average of the cross-validation errors close to zero (close is a relative term because it depends on the magnitudes of the data values), is the root mean square of the normalized errors close to one (normalized means divided by the kriging standard deviation). Look at a plot of the data locations coded by the cross validation errors and also a plot coded by the normalized cross validation errors. Look at the frequencies of normalized errors that are larger than 2.5 or less than -2.5, both frequencies should be small.
What evidence did you use to decide on an exponential model, try using a spherical model
- Badges
- Science method