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?


I couldn't answer your question using gstat package. However, you can also use geoR package to fit a variogram model to an empirical variogram and analyse fitted and residuals values. I give you a reproducible example below:

# Load libraries ----------------------------------------------------------


# Load data ---------------------------------------------------------------

#  Calcium content measured in soil data set

# Plot data

summary data plot

# Empirical omnidirectional variogram -------------------------------------

variogram <- variog(ca20)

# Variogram plot
plot(variogram, main = "Empirical variogram", pch = 19, col = "#0080FF")

enter image description here

# Fit model to empirical variogram ----------------------------------------

# Fit Spherical model
iniCovPars <- cbind("sigma" = seq(100, 200, length.out = 20), 
                    "phi" = seq(200, 600, length.out = 20))

# Using restricted maximum likelihood (REML) method
likfitVar.sph <- likfit(ca20, ini.cov.pars = iniCovPars, cov.model = "spherical", lik.method = "REML")

# Plot empirical variogram and adjusted model
plot(variogram, pch = 19, col = "#0080FF", ylim = c(0, 200), main = "Model fit to empirical variogram")
lines(likfitVar.sph, col = "red", lty = 1)
legend('topleft', lty = 1, legend = 'Spherical model', col = "red")

enter image description here

# Fitted
fitted <- fitted(likfitVar.sph)
fitted.df <- data.frame("east" = ca20$coords[,1], "north" = ca20$coords[,2], "fitted" = fitted)
coordinates(fitted.df) <- c("east","north")

# Residuals
res <- resid(likfitVar.sph)
res.df <- data.frame("east" = ca20$coords[,1], "north" = ca20$coords[,2], "fitted" = res)
coordinates(res.df) <- c("east","north")

# Original data
data.df <- data.frame("east" = ca20$coords[,1], "north" = ca20$coords[,2], "data" = ca20$data)
coordinates(data.df) <- c("east","north")

# Plots 
spplot(data.df, main = "Original data")

enter image description here

spplot(fitted.df, main = "Fitted values")

enter image description here

spplot(res.df, main = "Residuals values")

enter image description here

Note: check that you can compare differents models adjustments with AIC


[1] 1271.284
  • Hi, Guzmán, thanks for the detailed answer. However, this still is not the answer to my question. Fitting a variogram model to the empirical variogram basically looks like a simple (though non-linear) regression problem, in which the averaged semivariances depend on lag distances. In your third image you fitted a variogram model to your resulting empirical variogram. So for each lag-distance (the x-axis in the above plot) there is a fitted gamma value. If you have 15 lags, then there will be 15 semivaariances and hence 15 fitted values. – Asad Ali Jan 29 '17 at 11:51
  • @AsadAli Nice, I understand now what you were looking! I couldn't find an answer using the geoR package yet! I will update my answer when I find it. You can check your answer as the solution to your problem. – Guzmán Feb 2 '17 at 15:07

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

# 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.

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.