Software

On CRAN you can find the R package scpm (spatial smoothing with unknown change-points models). It allows you to perform 5 main tasks:

  1. (Residual) maximum likelihood estimation of covariance/semivariogram models for the spatial variation. For instance:
    • Mátern
    • Gneiting
    • Exponential
    • Gaussian, among others.
  2. (Residual) maximum likelihood estimation of traditional geostatistical linear models.
  3. Smoothing for geostatistical data or problems in 2D (lattices or scattered data).
    ATTENTION: for some scattered data or large datasets, estimation can be time demanding.
  4. Estimation of unknown change-points (1D), or contours of change (2D).
  5. Test to select between linear models and smoothing splines.

Download the package from https://cran.r-project.org/web/packages/scpm/index.html or within R console type install.packages("scpm"). An example code can be found below. The manual describes in detail the different commands and their theoretical aspects. If you have any comments, suggestions or questions you can contact me here.

		
#loading data landim1 originally from geoR package
data(landim1, package = "scpm")

library(scpm)

#converting data to class "sss"
d <- as.sss(landim1, coords = NULL, coords.col = 1:2, data.col = 3:4)

#Fitting spatial linear model with response A and covariate B
#Gneiting covariance function in the errors

#Null model
m0 <- scp(A ~ linear(~ B), data = d, model = "RMgneiting")

#Adding a bivariate cubic spline based on the coordinates
m1 <- scp(A ~ linear(~ B) + s2D(penalty = "cs"), data = d, model = "RMgneiting")

#Plotting observed and estimated field from each model
par(mfrow=c(2,2))
plot(m0, what = "obs", type = "persp", main = "Model null - y")
plot(m0, what = "fit", type = "persp", main = "Model null - fit")
plot(m1, what = "obs", type = "persp", main = "Model alternative - y")
plot(m1, what = "fit", type = "persp", main = "Model alternative - fit")

#Plotting the estimated semivariogram from each model
par(mfrow=c(1,2))
Variogram(m0, main = "Semivariogram - model null", ylim = c(0,0.7))
Variogram(m1, main = "Semivariogram - model alternative", ylim = c(0,0.7))

#Summary of the estimated coefficients
summary(m0)
summary(m1)

#Some information criteria
AIC(m0)
AIC(m1)
AICm(m0)
AICm(m1)
AICc(m0)
AICc(m1)
BIC(m0)
BIC(m1)