Regression with spatially varying coefficients in gstlearn
Xavier Freulon, Didier Renard, and the gstlearn Team
2025-04-23
Introduction
Regression models are build to explain a response variable over a region or domain. The regression coefficients are often assumed to be constant over the region of interest. However, in some applications, these coefficients may vary at the local level. The objective of this note is to present how this approach can be implemented using gstlearn.
Technical description
The variable of interest \(\{z(s), s \in \mathcal{D} \}\) is interpreted as the realization of some non stationary random function, \(Z(s)\) with \(Z(s) \in \mathbb{R}\) and \(\mathcal{D} \subset \mathbb{R}^d\) for \(d \leq 3\).
The model of kriging with external drifts
We consider first the case where the coefficients are constant but unknown: \[ Z(s) = \sum_{l = 0}^{L} a_l \, f^l(s) + Y(s). \] The coefficients of the spatial regression \(a_l\) are unknown and the residual \(Y\) is a centred stationary random function whose covariance function is \(C\).
The two first moments of \(Z\) are therefore:
\(m(s) = \mathbb{E}\{Z(s)\} = \sum_{l = 0}^{L} a_l \, f^l(s)\)
\(\text{Cov}(Z(s), Z(s')) = \mathbb{E}\{(Z(s)- m(s))(Z(s')- m(s'))\} = C(s - s')\)
The Best Linear Unbiased Predictor (BLUP) of \(Z\) at the target point \(\mathbf{o}\) is kriging with external drift usually considered in geostatistics
\[\hat{Z}(\mathbf{o}) = \sum_{\alpha \in A} \lambda^{\alpha}Z(\alpha)\]
with the kriging weights obtained minimizing the variance of the estimation error with the unbiasedness conditions.
The observations used to predict the variable \(Z\) at the target \(\mathbf{o}\) may be altered by some Gaussian additive noise (\(\epsilon \sim \mathcal{N}(0,1)\)):
\[\tilde{Z}(\alpha)= Z(\alpha) + \tau \, \varepsilon_{\alpha}\]
In this case, the kriging system can be written using matrix notation as
\[\left[ \begin{array}{cc} \mathbf{C}+ \tau^2 \mathbf{I} & \mathbf{F} \\ \mathbf{F}^T & \mathbf{0} \end{array} \right] \times \left[ \begin{array}{c} \Lambda \\ \mathbf{m} \end{array} \right] = \left[ \begin{array}{c} \mathbf{c} \\ \mathbf{d} \end{array} \right]\]
where the matrices and vectors used in the kriging system are:
the data covariance matrix, \(\mathbf{C}=[C(\alpha - \beta)]_{\alpha, \beta \in A}\)
the drift matrix at data locations, \(\mathbf{F}=[f^{l}(\alpha)]_{\alpha\in A, l \in 0:L}\) with \((f^0 = \mathbf{1})\),
the vector of the covariance between data and target, \(\mathbf{c}=[C(\alpha - \mathbf{o})]_{\alpha \in A}\),
the vector of the drifts at the target location, \(\mathbf{d}=[f^{l}(\mathbf{o})]_{l \in 0:L}\),
the vector of the kriging weights, \(\Lambda = [\lambda^{\alpha}]_{\alpha \in A}\),
the vector of the Lagrange coefficients , \(\mathbf{m}=[\mu_l]_{l \in 0:L}\).
The variance of the kriging error is given by: \[ \text{Var}(Z(\mathbf{o}) - \hat{Z}(\mathbf{o})) = C(0) - \Lambda^T \, \mathbf{c} - \mathbf{m}^T \, \mathbf{d} \]
With this model, different Best Linear Unbiased Estimators (BLUPs) can be considered to estimate different target variables. The first/left part of the kriging system is unchanged; only its second/right part must be modified to take into account the actual target variable. They are respectively,
- the estimation of the unobserved variable at target \(\mathbf{o}\), \(\hat{Z}(\mathbf{o})\), \[ \begin{array}{ccccc} \mathbf{c} & = & \left[\text{Cov}(\tilde{Z}(\alpha), Z(\mathbf{o}))\right]_{\alpha \in A} & = & [C(\alpha - \mathbf{o})]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [f^l(\mathbf{o})]_{l \in 0:L} \end{array}. \]
The variance of the variable to be estimated is \(\text{Var}(Z(\mathbf{o})) = C(0)\).
- the estimation of the residual at target \(\mathbf{o}\), \(\hat{Y}(\mathbf{o})\), \[ \begin{array}{ccccc} \mathbf{c} & = & \left[\text{Cov}(\tilde{Z}(\alpha), Y(\mathbf{o}))\right]_{\alpha \in A} & = & [C(\alpha - \mathbf{o})]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [0]_{l \in 0:L} \end{array}. \]
The variance of the residual to be estimated is \(\text{Var}(Z(\mathbf{o}) - m(\mathbf{o})) = C(0)\).
- the estimation of the drift at target \(\mathbf{o}\), \(\hat{m}(\mathbf{o}) = \hat{(\sum_{l\in 0:L} a_l\, f^l(0))}\), \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), \sum_{l\in 0:L} a_l\, f^l(\mathbf{o})\right]_{\alpha \in A} & = & [0]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [f^l(\mathbf{o})]_{l \in 0:L} \end{array} \]
The variance of the drift to be estimated is \(\text{Var}(m(\mathbf{o})) = 0\).
- the estimation of the unknown coefficients at target \(\mathbf{o}\) for \({l_0} \in 0:L\), \(\hat{a}_{l_0}]\), \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), a_{l_0}\right]_{\alpha \in A} & = & [0]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [\delta^l_{l_0}]_{l \in 0:L} \end{array} \]
The variance of the coefficient to be estimated is \(\text{Var}(a_{l_0}) = 0\).
In all cases, the variance of the estimation error of \(X\) is given by: \[ \text{Var}(X - \hat{X}) = Var(X) - \Lambda^T \, \mathbf{c} - \mathbf{m}^T \, \mathbf{d} \]
The linear model with spatially varying effects
A more complex model can be considered: some coefficients of the residual are varying (or not) across the domain of interest, the others are fixed. Latter on, the spatially varying coefficients are called the spatial effects, and the constant coefficients the fixed effects. The spatially effects are interpreted as the realization of some stationary random functions. The model for the variable of interest is therefore: \[ {Z}(\alpha) = \sum_{l = 0}^{P} A_l(\alpha) \, f^l(\alpha) + \sum_{l = P+1}^{L} a_l \, f^l(\alpha). \] The error \(\varepsilon \sim \mathcal{N}(0,\tau)\) may be added to the observations. This model includes the previous case:
we obtain the estimation of \(Z\) by kriging with external drift (KED) if \(P = 0\) and \(\tau = 0\) with \(L\) fixed effects,
we obtain the estimation of \(Z\) by ordinary kriging (OK) if \(L = 0\) and \(\tau = 0\), considering only the unknown mean.
More complex cases are obtained:
The estimation with spatially varying effects only if \(P = L\).
The mixed case with a residual, \(P\) spatially varying effects, and \(L-P\) fixed effect (\(0 < P < L\)).
If the varying coefficient \((A_l)_{l \in 0:P}\) are stationary random functions with unknown means \((a_l)_{l \in 0:P}\) and cross covariance functions \[\text{Cov}(A_l(s), A_p(s + \mathbf{h})) = C_{lp} (\mathbf{h})\] the model of the variable of interest is second order non stationary random function whose first two moments are:
\(m(s) = \mathbb{E}\{Z(s)\} = \sum_{l = 0}^{L} a_l \, f^l(s)\), and
\(C(s,s') = \text{Cov}(Z(s), Z(s')) = \sum_{l,l'\in 0:P} f^{l}(s) C_{l,l'}(s - s') f^{l'}(s')\).
Using this model, different BLUPs can be proposed, \[\hat{X} = \sum_{\alpha \in A} \lambda^{\alpha} \tilde{Z}(\alpha)\] The weights are defined by the kriging system but the matrix of the covariance between observations \(\mathbf{C}\) and the vector of the covariance between observations and the target \(\mathbf{c}\) are computed using non stationary covariance functions.
With this model, different Best Linear Unbiased Estimators (BLUPs) can be considered. The first part of the kriging system is unchanged and only its second part must be modified to take into account the actual target variable to be estimated. They are respectively:
- the estimation of the unobserved variable at target \(\mathbf{o}\), \(\hat{Z}(\mathbf{o})\) \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), Z(\mathbf{o}))\right]_{\alpha \in A} & = & [\sum_{l,l' \in 0:P} f^{l}(\alpha) C_{l,l'}(\alpha - \mathbf{o}) f^{l'}(0)]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [f^l(\mathbf{o})]_{l \in 0:L} \end{array} \]
The variance of the variable to be estimated is \(\text{Var}(Z(\mathbf{o})) = C(\mathbf{o},\mathbf{o})\).
- the estimation of the residual at target \(\mathbf{o}\), \(\hat{(Z(\mathbf{o}) - m(\mathbf{o}))}\) \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), Z(\mathbf{o}) - m(\mathbf{o}))\right]_{\alpha \in A} & = & [\sum_{l,l' \in 0:P} f^{l}(\alpha) C_{l,l'}(\alpha - \mathbf{o}) f^{l'}(0)]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [0]_{l \in 0:L} \end{array} \]
The variance of the variable to be estimated is \(\text{Var}(Z(\mathbf{o}) - m(\mathbf{o})) = C(\mathbf{o},\mathbf{o})\).
- the estimation of the drift at target \(\mathbf{o}\) \(\hat{m}(\mathbf{o}) = \hat{(\sum_{l\in 0:L} a_l\, f^l(0))}\) \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), \sum_{l\in 0:L} a_l\, f^l(\mathbf{o})\right]_{\alpha \in A} & = & [0]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [f^l(\mathbf{o})]_{l \in 0:L} \end{array} \]
The variance of the variable to be estimated is \(\text{Var}(m(\mathbf{o})) = 0\).
- the estimation of the spatial effects at target \(\mathbf{o}\) for \({l_0} \in 0:P\), \(\hat{A}_{l_0}(\mathbf{o})]\) \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), A_{l_0}(\mathbf{o})\right]_{\alpha \in A} & = & [\sum_{l \in 0:P} f^{l}(\alpha) C_{l,l'}(\alpha - \mathbf{o}) ]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [\delta^l_{l_0}]_{l \in 0:L} \end{array} \]
The variance of the variable to be estimated is \(\text{Var}(A_{l_0}(\mathbf{o})) = C_{l_0, l_0}(0)\).
- the estimation of the fixed* *effects** at target \(\mathbf{o}\) for \({l_0} \in (P+1):L\), \(\hat{a}_{l_0}]\), \[ \begin{array}{ccccc} \mathbf{c} & = &\left[\text{Cov}(\tilde{Z}(\alpha), a_{l_0}\right]_{\alpha \in A} & = & [0]_{\alpha \in A}\\ \mathbf{d} & = &\left[\sum_{\beta \in A} \lambda^{\beta} f^l(\beta) \right]_{l \in 0:L} & = & [\delta^l_{l_0}]_{l \in 0:L} \end{array} \]
The variance of the variable to be estimated is \(\text{Var}(a_{l_0}) = 0\).
As in previous section, the variance of the estimation error of the variable \(X\) is given by: \[ \text{Var}(X - \hat{X}) = Var(X) - \Lambda^T \, \mathbf{c} - \mathbf{m}^T \, \mathbf{d} \]
The equations of the kriging systems
In this section, we detail the estimation of the spatially varying effect \(A_{l_0}(\mathbf{o})\), for target \(\mathbf{o}\) and for \(l_0 \in 0:P\), using the Best Linear Unbiased Predicator (BLUP):
The linear predictor (LP) gives: \[ \hat{A}_{l_0}(\mathbf{o}) = \sum_{\alpha \in A} \lambda_{l_0}^{\alpha} \, \tilde{Z}(\alpha) \]
The unbiased condition (U) \(\mathbb{E}\{A_{l_0}(\mathbf{o}) - \hat{A}_{l_0}(\mathbf{o})\} = 0\) defines the \(L+1\) conditions (\(l \in 0:L\)): \[ \sum_{\alpha \in A} \lambda_{l_0}^{\alpha} \, f^l(\alpha) = \delta_{l_0}^l. \]
The optimal condition (B) implies that the variance of the estimation error be minimal under the unbiased conditions \[ \text{Var} \{A_{l_0}(\mathbf{o}) - \hat{A}_{l_0}(\mathbf{o})\} = C_{{l_0}{l_0}}(\mathbf{o}) - 2 \, \sum_{\alpha \in A} \lambda_{l_0}^{\alpha} c_{l_0}(\mathbf{o}, \alpha) + \sum_{\alpha, \beta \in A} \lambda_{l_0}^{\alpha} [C(\alpha, \beta) + \tau^2 \delta(\alpha - \beta)] \lambda_{l_0}^{\beta}, \] with the non stationary covariance functions \[ C(s, s') = \text{Cov} (Z(s), Z(s')) = \sum_{l,l' \in 0:P} f^{l}(s) \, C_{l,l'}(s - s') \, f^{l'}(s') \] and \[ c_l(s, s') = \text{Cov} (A_l(s), Z(s')) = \sum_{l' \in 0:P} C_{l,l'}(s - s') f^{l'}(s'). \]
The kriging system to solve is therefore: \[ \begin{array}{ccc} \sum_{\beta \in A} [C(\alpha, \beta) + \tau^2\delta(\alpha-\beta)] \, \lambda_{l_0}^{\beta} + \sum_{l \in 1:L} \mu_l \, f^l(\alpha) & = & c_{l_0}(\mathbf{o},\alpha) \\ \sum_{\beta \in A} f^l(\beta) \lambda_{l_0}^{\beta} & = & \delta_{l_0}^l \end{array} \]
In matrix form, the kriging system boils down to
\[\left[ \begin{array}{cc} \mathbf{C} + \tau^2 \mathbf{I} & \mathbf{F} \\ \mathbf{F}' & \mathbf{0} \end{array} \right] \times \left[ \begin{array}{c} \mathbf{\Lambda} \\ \mathbf{m} \end{array} \right]=\left[ \begin{array}{c} \mathbf{c} \\ \mathbf{d} \end{array} \right]\] where
the \([n, n]\) matrices defined by \(\mathbf{M} = [C(\alpha, \beta)]\) and \(\mathbf{I} = [\delta(\alpha - \beta)]\),
the \([n, L+1]\) matrix defined by \(\mathbf{F} = [f^l(\alpha)]\),
the \([n]\) vectors \(\mathbf{\Lambda} = [\lambda_{l_0}^{\alpha}]\) and \(\mathbf{c} = [c_{l_0}(\mathbf{o},\alpha)]\),
the \([L+1]\) vectors \(\mathbf{m} = [\mu_l]\) and \(\mathbf{d} = [\delta_{l_0}^l]\).
The number of observations is noted \(n = \#\{A\}\).
The number of linear predictors is noted \(L\).
The variance of the estimation error is: \[ \mathbf{Var} (A_{l_0}(\mathbf{o}) - \hat{A}_{l_0}(\mathbf{o})) = C_{l_0, l_0}(\mathbf{o},\mathbf{o}) - \sum_{\alpha \in A} \lambda_{l_0}^{\alpha} c_{l_0}(\mathbf{o},\alpha) - \mu_{l_0} = C_{l_0, l_0}(\mathbf{o},\mathbf{o}) - \Lambda^T \mathbf{c} - \mathbf{m}^T\mathbf{d} \]
The BLU predictor of the observed field \(Z\) requires to solve a similar kriging system where only the \([L+1]\) vector \(\mathbf{d} = [\delta_{l_0}^l \, f^{l}(\mathbf{o})]\) of the right part of the kriging system is modified.
The variance of the estimation error is: \[ \mathbf{Var} (Z(\mathbf{o}) - \hat{Z}(\mathbf{o})) = C(\mathbf{o},\mathbf{o}) - \sum_{\alpha \in A} \lambda_{l_0}^{\alpha} c_{l_0}(\mathbf{o},\alpha) -\sum_{l \in 0:L} \mu_l f^l(\mathbf{o}) = C(\mathbf{o},\mathbf{o}) - \Lambda^T \mathbf{c} - \mathbf{m}^T\mathbf{d} \]
Finally, the estimation of a fixed effect \(a_{l_0}\) (\(l_0 \in 0:L\)) requires to modify the right part of the kriging system with the \([n]\) vector \(\mathbf{c} = \mathbf{0}\) and the \([L+1]\) vector \(\mathbf{d} = [\delta_{l_0}^l]\).
The variance of this estimation error is therefore: \[ \mathbf{Var} (\hat{a}_{l_0}) = - \mu_{l_0} \]
Simulation of the reference case
We are defining some global environment variables within the gstlearn package.
set.seed(43243)
pal = rev(RColorBrewer::brewer.pal(6, "RdYlBu"))
opers = EStatOption_fromKeys(c("NUM", "MINI", "MAXI", "MEAN", "STDV"))
err = OptCustom_define("NotOptimSimpleCase", 1)
# list of estimations to be performed
ESTIM_KEY = c("obs", "res", "drift", "coeff")
estim_list = list()
estim_list[[1+length(estim_list)]] <- list(type = 1, idx = -1, name = "Z")
estim_list[[1+length(estim_list)]] <- list(type = 2, idx = -1, name = "RES")
estim_list[[1+length(estim_list)]] <- list(type = 3, idx = -1, name = "DRIFT")
estim_list[[1+length(estim_list)]] <- list(type = 4, idx = 0, name = "A")
estim_list[[1+length(estim_list)]] <- list(type = 4, idx = 1, name = "B")
suffix_key = c("estim", "stdev")
suffix_dsc = c("Estimation", "Std.")
flag_est = TRUE
flag_std = TRUE
flag_varz = FALSE
verbose = TRUEThe reference data set used to evaluate the different methods is simulated on the unit square \([0,1]^2\). The parameters are first simulated on a fine grid, then the observation are extracted from the grid. For each observation \(\alpha \in \mathcal{D}\), we get the observed value \(Z(\alpha)\) and the drifts \((f^l(\alpha))_{l \in 0:L}\). Note that the \(f^0\) is the constant value.
nx = c(50, 50)
dx = c(1,1) / nx
grd = DbGrid_create(x0 = c(0,0), dx = dx, nx = nx, angles = c(0, 0))
# Model of the spatially varying effects
tau = 0.1
m = c(1.25, 6.5)
rho = 0.7
sigma = c(2, 3)
c0 = diag(sigma^2)
c0[1,2] = prod(sigma)*rho
c0[2,1] = prod(sigma)*rho
model = Model_createFromParam(ECov_MATERN(), sills = c0, range = 0.5, param = 1)
err = model$setMeans(m)
# simulation of the spatially varying effects
err = simtub(dbin = NULL, dbout = grd, model = model, nbsimu = 1, nbtuba = 1000,
namconv = NamingConvention("Y"))
err = grd$setName(old_name = "Y.1", name = "A")
err = grd$setName(old_name = "Y.2", name = "B")
# simulation of the external drift Lognormal (m = 3, sigma = 0.25)
m_N <- 3
sigma_N <- 0.25
model_N = Model_createFromParam(ECov_CUBIC(), sill = 1.0, ranges = c(0.75, 0.75))
err = simtub(dbin = NULL, dbout = grd, model = model_N, nbsimu = 1, nbtuba = 1000,
namconv = NamingConvention("N"))
grd["N"] <- m_N * exp(sigma_N * grd["N"] - 1/2 * sigma_N^2)
# Computing the observed values
grd["DRIFT"] = grd["A"] + grd["B"]*grd["N"]
grd["RES"] = tau * rnorm(grd$getNSample())
grd["Z"] = grd["DRIFT"] + grd["RES"]
# names of the estimated variables
actual_var = c("Z", "RES", "DRIFT", "A", "B")
# a selection on the output grid (a disk)
grd["sel"] <- (sqrt((grd["x1"] - 0.5)^2 + (grd["x2"] - 0.5)^2) < 0.5)
# defining the observations
np <- 100
data = Db_createSamplingDb(grd, number = np, names = c("x1", "x2", "N", "Z"))
err = data$addSelectionRandom(prop = 0.8, seed = 13167, name = "sel")
# setting the roles
# on the grid
err = grd$setLocators("Z", ELoc_Z(), cleanSameLocator = TRUE)
err = grd$setLocators("N", ELoc_F(), cleanSameLocator = TRUE)
err = grd$setLocators("sel", ELoc_SEL(), cleanSameLocator = TRUE)
# on the data
err = data$setLocators("N", ELoc_F(), cleanSameLocator = TRUE)
err = data$setLocators("Z", ELoc_Z(), cleanSameLocator = TRUE)
err = data$setLocators("sel", ELoc_SEL(), cleanSameLocator = TRUE)
# the number of points in the moving neighborhood
n_neigh = 20
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| A | 1941 | -3.74 | 7.56 | 1.52 | 1.83 |
| B | 1941 | -2.88 | 17.48 | 7.46 | 3.73 |
| N | 1941 | 2.42 | 4.79 | 3.46 | 0.47 |
| Z | 1941 | -11.06 | 65.52 | 27.60 | 14.73 |
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| N | 80 | 2.36 | 4.73 | 3.29 | 0.52 |
| Z | 80 | -8.97 | 61.07 | 26.45 | 17.09 |
Kriging with external drift with a unique/moving neighbourhood
Initialisation
flag.unique.ked = FALSE
flag.unique = flag.unique.ked
if (flag.unique) {
neigh = NeighUnique()
prefix_section = "KED"
neigh_txt = "UNIQUE"
} else {
neigh = NeighMoving_create(nmaxi = n_neigh, nmini = n_neigh)
prefix_section = "KEDM"
neigh_txt = "MOVING"
}
prefix_ked = prefix_section
title_section = paste("using", neigh_txt, "neigh.")
modED = model$createReduce(0)
err = modED$addCovFromParam(ECov_NUGGET(), sill = tau)
err = modED$setDriftIRF(order = 0, nfex = 1)
define_estim <- function(estim, model, verbose = FALSE) {
# --------------------------------------------
# estimation of the observed variable
# --------------------------------------------
if (estim$type == 1) {
for (idx in 0:(model$getNCov()-1)) {
err = model$setCovFiltered(icov = idx, filter = FALSE)
}
for (idx in 0:(model$getNDrift()-1)) {
err = model$getDriftList()$setFiltered(i = idx, filter = FALSE)
}
}
# --------------------------------------------
# estimation of the residual
# --------------------------------------------
if (estim$type == 2) {
for (idx in 0:(model$getNCov()-1)) {
err = model$setCovFiltered(icov = idx, filter = FALSE)
}
for (idx in 0:(model$getNDrift()-1)) {
err = model$getDriftList()$setFiltered(i = idx, filter = TRUE)
}
}
# --------------------------------------------
# estimation of the observed drift
# --------------------------------------------
if (estim$type == 3) {
for (idx in 0:(model$getNCov()-1)) {
err = model$setCovFiltered(icov = idx, filter = TRUE)
}
for (idx in 0:(model$getNDrift()-1)) {
err = model$getDriftList()$setFiltered(i = idx, filter = FALSE)
}
}
# --------------------------------------------
# estimation of the coefficient
# --------------------------------------------
if (estim$type == 4) {
for (idx in 0:(model$getNCov()-1)) {
err = model$setCovFiltered(icov = idx, filter = TRUE)
}
for (idx in 0:(model$getNDrift()-1)) {
err = model$getDriftList()$setFiltered(i = idx, filter = !(idx == estim$idx))
}
}
# --------------------------------------------
if (verbose) {
model$display()
}
invisible()
} Estimation
# cleaning the target data base
err = grd$deleteColumns(names = paste0(prefix_section, ".*.", c("estim", "stdev"), ".*"))
for (iest in seq_along(estim_list)) {
if (verbose) {
print(paste0(">>> Estimation of ", estim_list[[iest]]$name, " using KED"))
}
# reset the model
err = define_estim(estim = estim_list[[iest]], model = modED, verbose = FALSE)
if (estim_list[[iest]]$type != 4) {
var_key = paste(prefix_section, ESTIM_KEY[estim_list[[iest]]$type], sep = ".")
} else {
var_key = paste(prefix_section, ESTIM_KEY[estim_list[[iest]]$type], estim_list[[iest]]$idx, sep = ".")
}
# computing the kriging
err = kriging(data, dbout = grd, model = modED, neigh = neigh, krigopt = KrigOpt(),
flag_est = flag_est, flag_std = flag_std,
namconv = NamingConvention(var_key))
# rescaling the drift component by the deterministic predictor to get the effect estimate
if ((estim_list[[iest]]$type == 4)&(estim_list[[iest]]$idx > 0)) {
if (flag_est) {
grd[paste(var_key, nm_obs, suffix_key[1], sep = ".")] =
grd[paste(var_key, nm_obs, suffix_key[1], sep = ".")] /
grd[nm_fac[estim_list[[iest]]$idx]]
}
if (flag_std) {
grd[paste(var_key, nm_obs, suffix_key[2], sep = ".")] =
grd[paste(var_key, nm_obs, suffix_key[2], sep = ".")] /
abs(grd[nm_fac[estim_list[[iest]]$idx]])
}
}
}## [1] ">>> Estimation of Z using KED"
## [1] ">>> Estimation of RES using KED"
## [1] ">>> Estimation of DRIFT using KED"
## [1] ">>> Estimation of A using KED"
## [1] ">>> Estimation of B using KED"Check
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| KEDM.obs.Z.estim | 1941 | -8.97 | 61.07 | 27.92 | 14.56 |
| KEDM.obs.Z.stdev | 1941 | 0.00 | 1.44 | 0.80 | 0.26 |
| KEDM.res.Z.estim | 1941 | -22.64 | 18.88 | -0.45 | 6.79 |
| KEDM.res.Z.stdev | 1941 | 1.06 | 1.84 | 1.49 | 0.13 |
| KEDM.drift.Z.estim | 1941 | 4.11 | 46.69 | 28.38 | 9.71 |
| KEDM.drift.Z.stdev | 1941 | 1.03 | 1.75 | 1.27 | 0.13 |
| KEDM.coeff.0.Z.estim | 1941 | -56.85 | 55.13 | -9.99 | 22.10 |
| KEDM.coeff.0.Z.stdev | 1941 | 2.82 | 5.95 | 4.28 | 0.77 |
| KEDM.coeff.1.Z.estim | 1941 | -3.20 | 22.25 | 11.15 | 4.85 |
| KEDM.coeff.1.Z.stdev | 1941 | 0.83 | 1.88 | 1.24 | 0.24 |
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| KEDM.diff.Z | 1941 | -13.26 | 10.35 | -0.32 | 3.27 |
| KEDM.diff.RES | 1941 | -18.71 | 22.37 | 0.46 | 6.79 |
| KEDM.diff.DRIFT | 1941 | -23.83 | 24.90 | -0.78 | 7.76 |
| KEDM.diff.A | 1941 | -50.01 | 58.26 | 11.51 | 20.98 |
| KEDM.diff.B | 1941 | -19.30 | 19.71 | -3.69 | 7.62 |
Kriging with spatially varying coefficients with a unique/moving neighbourhood
Initialisation
flag.unique.svc = TRUE
flag.unique = flag.unique.svc
if (flag.unique) {
neigh = NeighUnique()
prefix_section = "SVC"
neigh_txt = "UNIQUE"
} else {
neigh = NeighMoving_create(nmaxi = n_neigh, nmini = n_neigh)
prefix_section = "SVCM"
neigh_txt = "MOVING"
}
prefix_svc = prefix_section
title_section = paste("using", neigh_txt, "neigh.")
nfac = length(nm_fac)
model_svc = model$clone()
C0 = matrix(c(tau, 0, 0, 0), nrow = 2, ncol = 2)
err = model_svc$addCovFromParam(type = ECov_NUGGET(), sills = C0)
err = model_svc$setDriftIRF(order = 0, nfex = nfac)Estimation
# cleaning the target data base
err = grd$clearLocators(locatorType = ELoc_Z())
err = grd$deleteColumns(names = paste0(prefix_section, ".*.", c("estim", "stdev"), ".*"))
err = kriging_SVC(dbin = data, dbout = grd, model = model_svc, neigh = neigh,
flag_est = flag_est, flag_std = flag_std, flag_varz = flag_varz,
estim_list = estim_list,
prefix = prefix_section, verbose = verbose
)## [1] "SVC: number of spatial effects = 2"
## [1] "SVC: number of constant effects = 0"
## [1] "SVC: the spatially varying factors are not centered (2 universal conditions)"
## [1] "SVC: number of estimated variables = 10"
## [1] "SVC: variable #1 = obs.Z"
## [1] "SVC: variable #2 = res.Z"
## [1] "SVC: variable #3 = drift.Z"
## [1] "SVC: variable #4 = coeff.0.Z"
## [1] "SVC: variable #5 = coeff.1.Z"Checking
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| SVC.obs.Z.estim | 1941 | -8.97 | 61.07 | 27.87 | 14.45 |
| SVC.obs.Z.stdev | 1941 | 0.00 | 8.37 | 4.04 | 1.76 |
| SVC.res.Z.estim | 1941 | -31.67 | 31.15 | 0.77 | 13.29 |
| SVC.res.Z.stdev | 1941 | 3.71 | 10.97 | 6.39 | 1.52 |
| SVC.drift.Z.estim | 1941 | 17.33 | 39.69 | 27.10 | 4.48 |
| SVC.drift.Z.stdev | 1941 | 3.71 | 8.55 | 4.85 | 1.05 |
| SVC.coeff.0.Z.estim | 1941 | -10.44 | -0.77 | -5.53 | 1.85 |
| SVC.coeff.0.Z.stdev | 1941 | 10.37 | 11.32 | 10.89 | 0.19 |
| SVC.coeff.1.Z.estim | 1941 | 0.49 | 18.45 | 9.60 | 3.46 |
| SVC.coeff.1.Z.stdev | 1941 | 2.39 | 4.30 | 3.38 | 0.33 |
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| SVC.diff.Z | 1941 | -12.92 | 14.67 | -0.26 | 3.39 |
| SVC.diff.RES | 1941 | -31.24 | 31.40 | -0.77 | 13.29 |
| SVC.diff.DRIFT | 1941 | -32.20 | 37.27 | 0.50 | 13.66 |
| SVC.diff.A | 1941 | 3.76 | 9.91 | 7.05 | 1.03 |
| SVC.diff.B | 1941 | -5.33 | 1.72 | -2.15 | 0.92 |
KED vs. SVC
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| Z | 1941 | -11.06 | 65.52 | 27.60 | 14.73 |
| KEDM.obs.Z.estim | 1941 | -8.97 | 61.07 | 27.92 | 14.56 |
| SVC.obs.Z.estim | 1941 | -8.97 | 61.07 | 27.87 | 14.45 |
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| A | 1941 | -3.74 | 7.56 | 1.52 | 1.83 |
| KEDM.coeff.0.Z.estim | 1941 | -56.85 | 55.13 | -9.99 | 22.10 |
| SVC.coeff.0.Z.estim | 1941 | -10.44 | -0.77 | -5.53 | 1.85 |
| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| B | 1941 | -2.88 | 17.48 | 7.46 | 3.73 |
| KEDM.coeff.1.Z.estim | 1941 | -3.20 | 22.25 | 11.15 | 4.85 |
| SVC.coeff.1.Z.estim | 1941 | 0.49 | 18.45 | 9.60 | 3.46 |