1
$\begingroup$

I am having problems fitting a variogram model. I tried to change some parameters to estimate or fix them but I am still not achieving any improvement.

I remove trend of the data and use logarithms because of the asymmetry of the data.

There is some way to improve estimations and get a better fit for the variogram?

Data:

Meuse_zinc from the library(sp) and I choose zinc for y value (transform by logarithms because of the asymmetry of the data)

    library(sp)
    library(geoR)

    data(meuse) select<-c(1,2,6) meuse_zinc<-meuse[,select];
    logmeuse_zinc<-meuse_zinc 
    logmeuse_zinc$zinc<-log(meuse_zinc$zinc)
    loggeozinc<-as.geodata(logmeuse_zinc)

The likfit( ) function to estimate all parameters:

geodata<-loggeozinc
modelos<-c( "gaussian","exponential", "spherical","matern")
params<-c("Nugget", "Partial sill", "Range", "AIC", "Kappa")

modelo_lk<-function(geodata, modelo, kappa_value){
  ml1 <- likfit(geodata, 
trend=~poly(geodata$coords[,1], 2) + poly(geodata$coords[,2], 2), cov.model = modelo, 
                ini = c(0.86,0.19), fix.nugget =FALSE, nugget = 0.01, 
                fix.psiA = F, psiA = pi/2, fix.psiR = F, psiR = 3,
                fix.kappa = FALSE, kappa=kappa_value,
                lik.method = "REML" 
                )
  return(ml1)
}  

gaussian_lk<-modelo_lk(geodata, "gaussian", kappa_value = 0.5)
exponential_lk<-modelo_lk(geodata, "exponential", kappa_value = 0.5)
spherical_lk<-modelo_lk(geodata, "spherical", kappa_value = 0.5)
matern_lk<-modelo_lk(geodata, "matern", kappa_value = 0.3)

And I always get:

WARNING: estimated range is less than 1 tenth of the minimum distance between two points. Consider re-examine the model excluding spatial dependence 

And so when I plot them are practically equal:

maxdist <- variog(geodata, option = "cloud")$max.dist;
bin.1 <- variog(geodata, trend=~poly(geodata$coords[,1], 2)+poly(geodata$coords[,2],2),option = "bin", pairs.min=25, max.dist=maxdist/2,
                    estimator.type = "modulus");


plot(bin.1,  main="PARAMETRIC VARIOGRAMS", cex.main=1, pch=16);
lines(exponential_lk, lwd = 2, col="red", max.dist=maxdist/2);
lines(gaussian_lk, lwd = 2, col="blue", max.dist=maxdist/2);
lines(spherical_lk, lwd = 2, col="green3", max.dist=maxdist/2);
lines(matern_lk, lwd = 2, col="yellow", max.dist=maxdist/2);
legend(x="bottomright", inset=0.01, lty=c(1, 1), col=c("red", "blue", "green","yellow"),
       legend=c("Exponential", "Gaussian", "Spherical","Matern"), cex=1);

And also when I consider the envelope with confidence bands, there is a poor fit:

    env <- variog.model.env(geodata = geodata,coords=geodat$coords,obj.variog=bin.1, model.pars = exponential_lk);
plot(bin.1,main="CONFIDENCE BANDS FOR EMPIRICAL VARIOGRAM\nMODEL:Exponential", lwd=2,pch=16, 
     envelope = env);

That is how the empirical variogram looks: enter image description here

maxdist <- variog(loggeozinc, option = "cloud")$max.dist;
vario.b.robust <- variog(geodata=loggeozinc, option = "bin", pairs.min=30, max.dist=maxdist/2,estimator.type = "modulus")
plot(vario.b.robust, main = "EMPIRICAL VARIOGRAMS (MODULUS)\nBINNED",cex.main=1, cex.lab=1, cex=1, pch=16);

The variogram without trend (~poly(x, 2)+poly(y,2)) : enter image description here

bin <- variog(loggeozinc, trend=~poly(loggeozinc$coords[,1], 2)+poly(loggeozinc$coords[,2],2),
            option = "bin", 
            pairs.min=30, max.dist=maxdist/2,
            estimator.type = "modulus");
$\endgroup$
  • 3
    $\begingroup$ Variography is both (a) an art and (b) inherently exploratory. Please don't let the availability of automatic software fitting procedures persuade you otherwise. Thus, the very first thing to do in addressing this problem is to plot the experimental variogram. What does it look like? $\endgroup$ – whuber Oct 18 '18 at 15:56
  • $\begingroup$ I choose to do the robust variogram and seems there is spatial correlation, however the shape is quite strange, regarding the theorical variogram. I did explore the isotropy of the data and by doing rose diagram it seems there could be geometric anisotropy.That is the point now which model fits better this variogram. I am quite new in this area so maybe I did not make the better approach to solve that I appreciate your advice. I am going to try with ordinary least squares to estimate the parameters, but there could be a problem with likelihood or something I forget? $\endgroup$ – Nek Oct 19 '18 at 18:02
  • $\begingroup$ Why are you initializing the range to 0.19? The plot makes it clear that reasonable values must be orders of magnitude greater. Change 0.19 to 1000 or thereabouts and try again. Work on the spherical model. $\endgroup$ – whuber Oct 21 '18 at 21:23
  • $\begingroup$ Thank you, I was wrong thinking ini may be sill possible values (maximum, minimum), that helps me a lot. So ini is for (partial sill, range) possible reasonable values could be (0.3, 1000)? By doing that, the estimated range is 800 however estimated partial sill result is higher than 0.4, how could be posible? $\endgroup$ – Nek Oct 22 '18 at 13:38
  • $\begingroup$ The plot of the experimental variogram itself suggests the sill should be around 0.5 to 0.6. $\endgroup$ – whuber Oct 22 '18 at 15:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.