Multivariable prediction

When $s$ variables $Z_k (s)$, $k=1,...,m$ each follow a linear model $Z_k (s) = F_k \beta_k + e_k (s)$, and the $e_k (s)$ are correlated, then it makes sense to extend the weighted least squares model to allow multivariable prediction. Without loss of generality, assume $m=2$. When ${\bf z}(s)=(z_1 (s),z_2 (s))'$ and ${\bf B} = (\beta_1, \beta_2)'$ are substituted for $z(s)$ and $\beta$, and when

$$\begin{eqnarray*} {\bf f}(s_0) = \left[ \begin{array}{cc} f^1(s_0) & 0 \\ 0 & f... ...y}{cc} v_{11} & v_{12} \\ v_{21} & v_{22} \end{array} \right] \end{eqnarray*}$$

with $f^k(s_0)$ the $f(s_0)$ that corresponds to variable $k$, with

$$V_{21} = [{\rm Cov}(e_2 (s_i),e_1 (s_j))],$$

$$v_{21} = ({\rm Cov}(e_2 (s_1 ),e_1 (s_0 )),..., {\rm Cov}(e_2 (s_n ), e_1 (s_0 )))',$$

and 0 a conforming zero matrix or vector, are substituted for $f(s_0)$, $F$, $V$ and $v_0$, then the left-hand sides of both (A.2) and (A.4) yield the multivariable predictions: the left-hand side of (A.2) then becomes the prediction vector $\hat {\bf z} (s_0 ) = ( \hat{z}_1(s_0 ), \hat{z}_2 (s_0 ))'$, and the left-hand side of (A.4) becomes the ($2 \times 2$) matrix with prediction covariances.

The examples above assume that each variable $Z_k (s)$ has it's own unique parameter vector $\beta_k$. It is however possible to define common parameters for multiple variables (with merge). It is for instance possible to define a common mean (intercept) (see ``standardised ordinary cokriging'' in [8, p. 70], [12, p. 409-416], or [10, p. 323], or ``collocated cokriging'' [10, p. 326]), or to define a common regressor across different variables (cf. Analysis of covariance models, [3, p. 193 and excercise 9.1]). In any case, ${\bf f} (s_0 )$ and ${\bf F}$ loose their typical block structure. For instance, with $m=2$, and the trend of both $Z_1 (s)$ and $Z_2 (s)$ consist of one single intercept (merging the intercept for both variables) leads to an ${\bf F}$ matrix with only one column, filled with ones.

Uncorrelated least squares prediction

When the residuals are not correlated and have unequal variance, i.e. under the model

$$Z(s) = F \beta + e(s), \hspace{.5cm}{\rm E}(e(s))= 0, \hspace{.5cm}{\rm Cov}(e(s)) = \sigma^2 D$$

with $D$ a known diagonal matrix, the uncorrelated least squares estimate for $\beta$ is obtained by

$$\hat{\beta}_* = (F'D^{-1} F)^{-1} F'D^{-1}z(s)$$

having variance

$${\rm Var}(\beta - \hat{\beta}_{*} ) = (F'D^{-1} F)^{-1} \sigma^2$$

where $\sigma^2$ is estimated by

$$s^2 = z(s)'(I-F D^{-1} (F'D^{-1} F)^{-1} F'D^{-1})z(s)/(n-p)$$

Least squares predictions at location $s_0$ are obtained by

$$\hat Z (s_0 ) = f(s_0 ) \hat{\beta}_{*}$$ (A5)

having variance

$${\rm Var}(Z (s_0 ) - \hat Z (s_0 )) = (1 + f(s_0 )(F'D^{-1} F)^{-1} f(s_0 )') \sigma^{2}$$ (A6)

or, when block prediction is involved

$${\rm Var}(Z (B_0 ) - \hat Z (B_0 )) = f(s_0 )(F'D^{-1} F)^{-1} f(s_0 )' \sigma^2$$ (A7)

When variograms are specified, trend prediction (using method: trend;) involves the calculation of $f(s_0 ) \hat \beta$, having variance $f (s_0 )(F'V^{-1} F)^{-1} f(s_0 )'$. When no variograms are specified, trend prediction involves the evaluation of (A.5) and for variances either (A.6) or, for blocks, (A.7). $D$ is by default the identity matrix $I$, in which case ordinary least squares (OLS) is used, or else it has elements specified by the V column of the data( id) file concerned. Specifying V also results in weighted uncorrelated least squares estimation of residuals in case of variogram modelling.

When common parameters are defined with merge, and no variograms are specified, then estimation and prediction under the OLS and WLS models will be done, in which case the constant variance is assumed for the joint (multivariable) residual. In this case, known ratios in variances between different variables can still be set using the V field of each data variable.

Kriging data with known measurement errors

Known, constant measurement error can be defined in the variogram model. Suppose the variable $y$ has an apparent nugget effect of 1 and a spatially correlated part of 1 exp(10), than the variogram model can be written as 1 Nug() + 1 Exp(10). This would yield the `standard' kriging predictions, i.e. exact interpolation. If it is known that 75% of the nugget variance constitues of measurement error, and predictions are required for the measurement error free part of the variable, then the variogram of $y$ can be defined as:

variogram(y): 0.75 Err() + 0.25 Nug() + 1 Exp(10)

see for more information [5].

Known, varying measurement errors can be defined (for data in column files) by specifying the V field. When variograms are specified and the goal is prediction or trend prediction, then the covariances ${\rm Cov}(e(s_i ), e(s_i ))$ are taken to be $c - \gamma(h) + \sigma^2_{\epsilon} (s_i )$, with $\sigma^2_\epsilon(s_i)$ the value of the V-field of record $i$, thus interpreting the variance field (V) as a known, location-specific measurement error [7,20,5]. Otherwise formulated: in this case $V+D$ is used as covariance matrix, instead of $V$ only. Putting a constant in the V field should yield the same results as specifying this value in the Err term of the variogram.