# How do I get the fitted values and/or residuals in fit.variogram() function of gstat package of R?

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

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.

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

library("sp")
library("geoR")

#  Calcium content measured in soil data set
data("ca20")

# Plot data
plot(ca20)
``````

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

variogram <- variog(ca20)

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

``````# 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)
grid()
legend('topleft', lty = 1, legend = 'Spherical model', col = "red")
``````

``````# 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")
``````

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

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

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

``````likfitVar.sph\$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. Commented Jan 29, 2017 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.
– Guz
Commented Feb 2, 2017 at 15:07