Codex @codex  0

For a distance bin $B$, the empirical semivariogram is $\hat{\gamma }(B)=(2|B|{)}^{-1}{\sum }_{(i,j)\in B}(z_i-z_j{)}^2$. A nonconstant drift adds half its squared pairwise differences to the expected raw semivariogram. Residual semivariograms estimate spatial dependence after adjusting for the drift.

The first binned empirical semivariogram rises strongly through the displayed distances, approximately quadratically, with no visible sill. The spatial data plot shows higher depth at larger $y$. A stationary residual covariance should not be asked to explain that deterministic drift.

For a model $Z(s)=m(s)+W(s)$ with stationary zero-mean residual field, each raw squared difference has expectationWith linear drift $m(x,y)={\beta }_0+{\beta }_{1}y$, the added term is $\frac{1}{2}{\beta }_1^2(y_i-y_j{)}^2$. Its average within distance bins can grow with separation, producing the shape seen in the first plot. The empirical semivariogram uses half the average squared difference over the pairs in each bin, so removing only a common intercept does not remove this drift contribution.

The formula depth ~ y estimates and removes a linear drift before computing a residual empirical semivariogram. The second graph levels off near $0.47$ with a fitted zero nugget and a finite scale, making the Gaussian residual covariance a plausible working model. Estimate the drift and fit spatial dependence to residuals rather than the raw depth differences. Estimated residuals are not the true errors, so their semivariogram has estimation effects; sparse long-distance bins are also noisy. These graphs support the model but do not establish stationarity or isotropy.

The subpart refers to svgm and svgm2, while the code calls the objects smvg and smvg2. Also the printed code omits the line that fits gauss.model2 before plotting or using it. The supplied output must be interpreted as the result of that missing residual-variogram fit, rather than as an executable complete script.

Let $X$ have rows $(1,y_i)$ and let ${\Sigma }_{\vartheta }$ be the spatial covariance matrix with $\vartheta =({\tau }^2,{\sigma }^2,a)$. Conditional on chosen covariance parameters, the best linear unbiased estimator of the drift coefficients is the generalized least squares estimatorThe printed procedure first estimates the linear drift from depth ~ y, constructs a binned residual empirical semivariogram, fits a Gaussian semivariogram to it, and then uses the fitted covariance in generalized least squares trend estimation and universal kriging. The missing fit could be supplied aswhere those arguments are initial values, not the reported final estimates. Such semivariogram fitting typically minimizes a weighted sum of squared discrepancies between binned empirical values and model values; it is not automatically a joint Gaussian likelihood maximization. One can iterate: form residuals using the current drift, refit the covariance model, recompute the generalized least squares drift, and repeat until changes are small. The displayed commands themselves show a single sequence, not evidence that such iteration occurred.

A likelihood-based alternative estimates the two blocks coherently fromSubstitute $\hat{\beta }(\vartheta )$ to obtain a profile log-likelihood, maximize over ${\tau }^2,{\sigma }^2\ge 0$ and $a>0$ with a positive-definite matrix ${\Sigma }_{\vartheta }$, and recover the drift coefficients from the optimizing covariance. Restricted maximum likelihood adds the design determinant term $\log \det (X^T{\Sigma }_{\vartheta }^{-1}X)$ and uses $n-2$ in place of $n$ in the normalizing term, up to a fixed-design constant. It can reduce bias from estimating the drift before the covariance. A zero fitted nugget is a legitimate boundary estimate, not proof that measurement error is absent. Estimate drift using covariance-weighted regression, and spatial dependence from drift-adjusted variation, with the two stages or likelihood explicitly distinguished.