sm.variogram: Confidence intervals and tests based on smoothing an empirical variogram.

Description

This function constructs an empirical variogram, using the robust form of construction based on square-root absolute value differences of the data. Flexible regression is used to assess a variety of questions about the structure of the data used to construct the variogram, including independence, isotropy and stationarity. Confidence bands for the underlying variogram, and reference bands for the independence, isotropy and stationarity models, can also be constructed under the assumption that the errors in the data are approximately normally distributed.

Usage

sm.variogram(x, y, h, df.se = "automatic", max.dist = NA, original.scale = TRUE,  varmat = FALSE, ...)

Arguments

a vector or two-column matrix of spatial location values.

a vector of responses observed at the spatial locations.

a smoothing parameter to be used on the distance scale. A normal kernel function is used and h is its standard deviation. However, if this argument is omitted h will be selected by an approximate degrees of freedom criterion, controlled by the df parameter. See sm.options for details.

df.se

the degrees of freedom used when smoothing the empirical variogram to estimate standard errors. The default value of "automatic" selects the degrees of smoothing described in the Bowman and Crujeiras (2013) reference below.

max.dist

this can be used to constrain the distances used in constructing the variogram. The default is to use all distances.

original.scale

a logical value which determines whether the plots are constructed on the original variogram scale (the default) or on the square-root absolute value scale on which the calculations are performed.

varmat

a logical value which determines whether the variance matrix of the estimated variogram is returned.

...

other optional parameters are passed to the sm.options function, through a mechanism which limits their effect only to this call of the function. An important parameter here is model which, for sm.variogram, can be set to "none", "independent", "isotropic" or "stationary". Other relevant parameters are add, eval.points, ngrid, se, xlab, ylab, xlim, ylim, lty; see the documentation of sm.options for their description. See the details section below for a discussion of the display and se parameters in this setting.

Value

sqrtdiff, distance: the raw differences and distances
sqrtdiff.mean, distance.mean: the binned differences and distances
weights: the frequencies of the bins
estimate: the values of the estimate at the evaluation points
eval.points: the evaluation points
h: the value of the smoothing parameter used
ibin: an indicator of the bin in which the distance between each pair of observations was placed
ipair: the indices of the original observations used to construct each pair
p: the p-value of the test
se: the standard errors of the binned values (if the argument se was set to TRUE)
se.band: when an independence model is examined, this gives the standard error of the difference between the smooth estimate and the mean of all the data points, if a reference band has been requested
V: the variance matrix of the binned variogram. When model is set to "isotropic" or "stationary", the variance matrix is computed under those assumptions.
sdiff: the standardised difference between the estiamte of the variogram and the reference model, evaluated at eval.points
levels: the levels of standarised difference at which contours are drawn in the case of model = "isotropy".

Side Effects

a plot on the current graphical device is produced, unless the option display="none" is set.

Details

The reference below describes the statistical methods used in the function. Note that, apart from the simple case of the indpendence model, the calculations required are extensive and so the function can be slow.

The display argument has a special meaning for this function. Its default value is "binned", which plots the binned version of the empirical variogram. As usual, the value "none" will suppress the graphical display. Any other value will lead to a plot of the individual differences between all observations. This will lead to a very large number of plotted points, unless the dataset is small.

References

Diblasi, A. and Bowman, A.W. (2001). On the use of the variogram for checking independence in a Gaussian spatial process. Biometrics, 57, 211-218. Bowman, A.W. and Crujeiras, R.M. (2013). Inference for variograms. Computational Statistics and Data Analysis, 66, 19-31.

Examples

Run this code

with(coalash, {
   Position <- cbind(East, North)
   sm.options(list(df = 6, se = TRUE))

   par(mfrow=c(2,2))
   sm.variogram(Position, Percent, original.scale = FALSE, se = FALSE)
   sm.variogram(Position, Percent, original.scale = FALSE)
   sm.variogram(Position, Percent, original.scale = FALSE, model = "independent")
   sm.variogram(East,     Percent, original.scale = FALSE, model = "independent")
   par(mfrow=c(1,1))
})

## Not run: 
# # Comparison of Co in March and September
#    
# with(mosses, {
# 	
#    nbins <- 12
#    vgm.m <- sm.variogram(loc.m, Co.m, nbins = nbins, original.scale = TRUE,
#                         ylim = c(0, 1.5))
#    vgm.s <- sm.variogram(loc.s, Co.s, nbins = nbins, original.scale = TRUE,
#                         add = TRUE, col.points = "blue")
#                         
#    trns <- function(x) (x / 0.977741)^4
#    del <- 1000
#    plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "b",
#          ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
#    points(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean), type = "b",
#          col = "blue", pch = 2, lty = 2)
# 
#    plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "b",
#          ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
#    points(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean), type = "b",
#          col = "blue", pch = 2, lty = 2)
#    segments(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean - 2 * vgm.m$se),
#          vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean + 2 * vgm.m$se))
#    segments(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean - 2 * vgm.s$se),
#          vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean + 2 * vgm.s$se),
#          col = "blue", lty = 2)
# 
#    mn <- (vgm.m$sqrtdiff.mean + vgm.s$sqrtdiff.mean) / 2
#    se <- sqrt(vgm.m$se^2 + vgm.s$se^2)
#    plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "n",
#         ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
#    polygon(c(vgm.m$distance.mean, rev(vgm.m$distance.mean)),
#         c(trns(mn - se), rev(trns(mn + se))),
#         border = NA, col = "lightblue")  
#    points(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean))
#    points(vgm.s$distance.mean, trns(vgm.s$sqrtdiff.mean), col = "blue", pch = 2)
# 
#    vgm1 <- sm.variogram(loc.m, Co.m, nbins = nbins, varmat = TRUE, 
#                         display = "none")
#    vgm2 <- sm.variogram(loc.s, Co.s, nbins = nbins, varmat = TRUE,
#                         display = "none")
# 
#    nbin  <- length(vgm1$distance.mean)
#    vdiff <- vgm1$sqrtdiff.mean - vgm2$sqrtdiff.mean
#    tstat <- c(vdiff %*% solve(vgm1$V + vgm2$V) %*% vdiff)
#    pval  <- 1 - pchisq(tstat, nbin)
#    print(pval)
# })
# 
# # Assessing isotropy for Hg in March
# 
# with(mosses, {
#    sm.variogram(loc.m, Hg.m, model = "isotropic")
# })
# 
# # Assessing stationarity for Hg in September
# 
# with(mosses, {
#    vgm.sty <- sm.variogram(loc.s, Hg.s, model = "stationary")
#    i <- 1
#    image(vgm.sty$eval.points[[1]], vgm.sty$eval.points[[2]], vgm.sty$estimate[ , , i],
#          col = topo.colors(20))
#    contour(vgm.sty$eval.points[[1]], vgm.sty$eval.points[[2]], vgm.sty$sdiff[ , , i],
#          col = "red", add = TRUE)
# })
# 
# ## End(Not run)

Run the code above in your browser using DataLab