# How can I understand these variograms?

Cross Validated Asked by Zheyuan Li on January 7, 2022

Using grf function from R package geoR, I simulated 6 replicates (each with 1000 samples) of a Gaussian random field on [0, 1] x [0, 1], with zero mean, zero nugget, and an exponential spatial covariance $$gamma(h) = expleft(frac{h}{phi}right)sigma^2,$$ with $sigma^2 = 1$ and $phi = 0.5$. In the following figure, the blue solid curve gives the semi-variogram of this true model and the blue dashed line gives the sill (which is $sigma^2 = 1$). Then for each replicate, I

1. compute an empirical semi-variogram with function variog (black solid dots);
2. compute a monte-carlo envelope for the empirical semi-variogram (gray shaded polygon);
3. fit an exponential spatial covariance model (black solid curve) to the empirical semi-variogram for $hat{sigma}^2$ and $hat{phi}$ (printed in the tile).

However, I find it difficult to understand the results / behavior as illustrated by the figure.

1. I simulated data from a true model and fit that true model, why does the fitted semi-variogram differ from the true one so much? Or, why does $(hat{sigma}^2, hat{phi})$ differ from $(sigma^2, phi)$ so much? On some replicate the difference is drastic.
2. How can the empirical semi-variogram (black solid dots), as well as the fitted variogram (black solid curve), breach sample variance (black dashed line)? Isn’t sample variance the sill?
3. Why does the black dashed line so much away from the blue dashed line? (the answer is maybe as same as that to question 1)

Is a spatial process expected to behave like this? If I examine an AR(1) process by simulating a time series and fitting it using Yule-Walker method, I can get back to the truth, i.e., true autocorrelation coefficient and variance. This is not possible for a spatial process?

Reproducible R code

set.seed(0)

phi <- 0.5
sigmasq <- 1

## unconditional simulation of Gaussian random field with exponential covariance
EXP_GRF <- grf(1000, cov.pars = c(sigmasq, phi), cov.model = "exponential",
nsim = 6)

## estimation of semi-variogram
semi_variog <- variog(EXP_GRF)

semi_variog_i <- semi_variog
semi_variog_i$v <- NULL par(mfrow = c(2, 3), mar = c(2.5, 2.5, 1, 1)) for (i in 1:6) { ## the i-th empirocal semi-variogram semi_variog_i$v <- semi_variog$v[, i] ## mc envelops env <- variog.mc.env(coords = EXP_GRF$coords, data = EXP_GRF$data[, i], obj.variog = semi_variog_i, nsim = 99, messages = FALSE) ## plotting range MAX <- max(env$v.lower, env$v.upper, semi_variog_i$v, sigmasq)
with(semi_variog_i, plot(u, v, ylim = c(0, MAX), type = "n"))

## empirical semi-variogram
polygon(c(semi_variog_i$u, rev(semi_variog_i$u)),
c(env$v.lower, rev(env$v.upper)), col = 8, border = NA)
with(semi_variog_i, points(u, v, pch = 19))

## variance of the data
var_data <- semi_variog$var.mark[i] abline(h = var_data, lty = 2) ## semi-variogram of the true model lines.variomodel(EXP_GRF, col = 4, lwd = 2) abline(h = sigmasq, lwd = 2, col = 4, lty = 2) ## fit an variogram for parameter estimation model <- variofit(semi_variog_i, c(var_data, phi), "exponential", fix.nugget = TRUE, nugget = 0, message = FALSE) ## fitted variogram lines.variomodel(model) ## add estimated sigmasq and phi as title title(sprintf("sigmasq = %.2f, phi = %.2f", model[[2]][1], model[[2]][2])) }  ## One Answer I didn't have a chance to look carefully into the code, but I can already guess what you are experiencing. It may have to due with the fact that you are simulating random fields on finite domains (in a computer). Try to increase the domain to a larger box [0,10]x[0,10] or try to decrease the range of the variogram to something that is at least smaller than half of your domain size, say 0.1. Then plot the variograms up to a lag that is at most half of your domain. Everything that is separated by a lag that is equal or larger to half of your domain is not useful in empirical variogram calculations (too few samples, and hence spurious oscillations). Geostatistics theory is more general than AR models, and there are theorems related to the issue you are experiencing, but I have no time at the moment to write down a more detailed answer. Answered by juliohm on January 7, 2022 ## Add your own answers! ## Related Questions ### R: When do we use mean or median for the y axis in ggplot2 when doing analysis on property prices? 0 Asked on January 28, 2021 by chua-s-yang ### COCO evaluation – Negative values on AP and AR 0 Asked on January 28, 2021 by visionenthusiast ### How to make the regressor of LASSO consistent? 0 Asked on January 28, 2021 by zqq ### Suggestions for identifying the most “important” image labels 1 Asked on January 28, 2021 by nlapidot ### Any ideas on how to segment a 2D vector field? 0 Asked on January 28, 2021 by tricostume ### Binomial logistic regression for multiclass problems 1 Asked on January 27, 2021 by mathews24 ### How is confidence defined in Expected Calibration Error? 0 Asked on January 26, 2021 by thecity2 ### Why does the McNemar’s test use$chi^{2}\$ and not the normal distribution?

2  Asked on January 26, 2021

### What algorithm can you use if you want clusters but only are interested in one group?

0  Asked on January 26, 2021 by bonesones

### Can I use an unknown number of variables to model my time-series?

0  Asked on January 26, 2021 by kplauritzen

### Variance of a stationary AR(2) model

2  Asked on January 26, 2021 by user369210

### Avoiding adjustments for time-varying controls in difference-in-differences (DID)?

0  Asked on January 26, 2021

### Removing the effect from structural breaks

1  Asked on January 25, 2021 by kiril-e-proykov

### Recommender System – Predict ratings with Random Forest Regressor or Classifier?

0  Asked on January 24, 2021 by oja-niva

### Nonparametric assessment of multiple predictors

0  Asked on January 24, 2021 by mephisto73

### Calculating measurement variance to achieve desired accuracy in estimation

0  Asked on January 23, 2021 by valjean

### Can large # of epochs or smaller batchsize compensate for smaller data size in training lstms

1  Asked on January 23, 2021 by tjt

### Probability that number of heads exceeds sum of die rolls

5  Asked on January 23, 2021 by user239903

### Combining Sub-Samples for Factor Analysis?

0  Asked on January 22, 2021

### Need to create a model to identify patterns in user details

0  Asked on January 21, 2021 by pooza