[GIS] How to get the fitted values and/or residuals in fit.variogram() function of gstat package of R

gstatrvariogram

I am trying to do some diagnostics after fitting a variogram model to my empirical variogram using fit.variogram() function in gstat package of R. I need to get a hold of residuals and the fitted values so that I can do normality and GOF tests. However, the fit.variogram() only provides the parameters' values (e.g. psill, nugget and rannge). How do I get residuals and/or the fitted values after fitting a variogram model?

Best Answer

I found a very simple solution for obtaing the fitted gamma values within gstat package. It's the variogramLine() function. A simple code is attached here.

# Empirical variogram
ev = variogram("pH", data=data,....)
fv = fit.variogram(q, "Sph", ....)
fitted=variogramLine(v.fit, maxdist=max(q$dist), dist_vector=q$dist)

fitted # see what are the values of fitted 

residual = fitted$gamma - q$gamma
sqrt(mean(residual^2))

# A plot will tell you what is happening..
plot(ev$dist, ev$gamma, type="b", cex=1, lty=1, ylim=c(0, 150000), main= "(a)", xlab="distance (meters)", ylab="semivariance",pch=16, col="darkgrey")
lines(fitted, type="b")

You get a plot like following.

Empirical variogram fitted with a spherical model at default 15 lags in gstat package

You see for each empirical gamma there is a fitted gamma value. This is what I wanted. I am actually after finding a the best fit model using some simple statistics.