CODEWITHSUNDEEP
rgeor
ToNoY
ToNoY

Reputation: 1378

fitting semivariogram using R

I've fitted my semivariogram using "geoR" package as follows:

#loading geoR
library(geoR)
#calculating the semivariogram
bin1 <- variog(geodata)
#plotting
plot(bin1)
#fitting the model
lines.variomodel(cov.model = "exp", cov.pars = c(0.571,0.2527), nug = 0.13, max.dist=1)

Is there any way to know the goodness-of-fit of the exponential model (with ease)??

Thanks in advance ...

Upvotes: 1

Views: 4780

Answers (1)

mnel
mnel

Reputation: 115390

You aren't really fitting the model here, you can use variofit to fit a model to the empirical semivariogram.

eg

library(geoR)
vario100 <- variog(s100, max.dist=1)
ini.vals <- expand.grid(seq(0,1,l=5), seq(0,1,l=5))
ols <- variofit(vario100, ini=ini.vals, fix.nug=TRUE, wei="equal")
## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: matern with fixed kappa = 0.5 (exponential)
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  1.1070  0.4006 
## Practical Range with cor=0.05 for asymptotic range: 1.200177
## 
## variofit: minimised sum of squares = 0.1025

You could calculate the appropriate weighted sum of squares your self at each of the points on the variogram using cov.spatial.

I don't think this is a good idea.

Instead, you could use loglik.GRF to calculate the likelihood associated with a given model and all the data

# you can pass an existing model
loglik.GRF(s100, obj.model = ols)
## [1] -87.32958

Or just the parameters

 loglik.GRF(s100, cov.pars = c(1.5,0.6), nugget = 0.01)

If you really want the weighted least squares value, you can use fit.variogram from the gstat package, without actually fitting anything.

One good thing is that gstat has some helper functions to convert from geoR models to gstat versions. This is really useful as gstat is much faster and more efficient for prediction.

eg

 library(gstat)
 # convert geodata to  data.frame object
 s100df <- as.data.frame(s100)
 # remove geodata.frame class that causes problems
 class(s100df) <- 'data.frame'

 # create gstat version of variogram
 s100v <- variogram(data~1, ~X1+X2, s100df)

 # convert a variomodel to vgm object

 foo <- as.vgm.variomodel(list(cov.model = 'exponential', kappa = 0.5,
                          cov.pars = c(1.5,0.6), nugget = 0.2))

 # get the weighted least squares value
 # calling fit.variogram without fitting any thing
 fittedfoo <- fit.variogram(s100v, foo, fit.sills = FALSE, fit.ranges = FALSE)
 # the weighted sum of squares is
 attr(fittedfoo, 'SSErr')
 ## [1] 0.6911813

Upvotes: 2

Related Questions