Decorrelation and Imputation Methods for Multivariate Modeling

Abstract

In most mining projects, multivariate modeling of regionalized variables has a critical impact on the final model due to complex multivariate relationships between correlated variables. In geostatistical modeling, multivariate transformations are commonly employed to model complex data relationships. This decorrelates or makes the variables independent, which enables the generation of independent models for each variable while maintaining the ability to restore multivariate relationships through a back-transformation. There are a myriad of transformation methods, however, this chapter discusses the most applied methods in geostatistical procedures. These include principal component analysis (PCA), minimum/maximum autocorrelation factors (MAF), stepwise conditional transform (SCT), and projection pursuit multivariate transform (PPMT). All these transforms require equally sampled data. In the case of unequal sampling, it is common practice to either exclude the incomplete samples or impute the missing values. Data imputation is recommended in many scientific fields as removing incomplete samples usually removes valuable information from modeling workflows. Three common imputation methods are discussed in this chapter: single imputation (SI), maximum likelihood estimation (MLE), and multiple imputation (MI). Bayesian updating (BU) is also discussed as an adaptation of MI to geostatistical analysis. MI methods are preferred in geostatistical analysis because they reproduce the variability of variables and reflect the uncertainty of missing values.

1. Introduction

Geostatistical modeling problems often consider more than one variable of interest. In the geosciences, these are obtained through log, core, and remote sensing measurements made in a subsurface domain of interest. Geostatistical modeling requires spatial characterization of multiple attributes to ensure the final model correctly reproduces the relations between all variables. Traditional geostatistical approaches that account for multivariate heterotopic data are based on a cokriging framework [1], which makes use of the linear model of coregionalization (LMC). To utilize LMC, all direct and cross-covariances between experimental variables must be modeled simultaneously. This can be a tedious task because some datasets may contain many variables, but automatic fitting software usually fits LMCs automatically.

To avoid the complexities of cokriging and LMC, simplified methodologies - such as Markov models [2, 3] and decorrelation methods - are used in spatial modeling workflows. Markov models impose the presence of an exhaustive secondary attribute that is used as auxiliary information to model a variable of interest [4]. Although Markov models are simpler methods for multivariate modeling, they are limited in their ability to reproduce the targeted spatial and multivariate relationships. Alternatively, decorrelation methods—such as principal components analysis (PCA) [5, 6], minimum/maximum autocorrelation factors (MAF) [7, 8], stepwise conditional transform (SCT) [9, 10], and projection pursuit multivariate transform (PPMT) [11, 12]—do not require the use of an LMC and cokriging. Decorrelation methods transform correlated variables into a set of uncorrelated/independent components or factors. Each factor can then be modeled independently, and the correlations between the attributes are restored through the back transform.

A major limitation of such methodologies is the need for homotopic sampling, which requires the same number of samples to be collected for each variable at all sampling locations [13]; in most practical applications it is challenging to adhere to homotopic sampling. This difficulty, rather than homotopic sampling itself, can lead to distortions in the known distributions of the variables, which could negatively impact the final model. Thus, heterotopic sampling configurations must be preprocessed before these multivariate techniques can be applied [14].

During data preprocessing when considering decorrelation workflows, samples with missing values in multivariate data are either excluded or imputed. Excluding samples with missing values is rarely considered because it reduces the amount of information available to build the model [15]; exclusion is only appropriate when there are few samples with missing values and the modeler decides that their removal will not have a negative impact. Utilizing data imputation maintains the initial dataset size; however, the approach must account for data variability, local conditional statistics, and spatial correlations in geostatistical analysis to ensure the final model is in alignment with the observed phenomena. It should be noted that conventional imputation methods are not specifically optimized for geostatistical analysis. Thus, imputation methods that consider spatial data, such as Bayesian updating (BU), are preferred in geological data analysis [16, 17].

This chapter introduces the fundamentals of multivariate transforms and missing data theory. The following section begins with an overview of multivariate transform theory and discusses common decorrelation methods such as PCA, MAF, SCT, and PPMT. After outlining these multivariate transforms, an overview of missing data theory is provided. The chapter concludes with a review of conventional missing data imputation methods—single imputation (SI), maximum likelihood estimation (MLE), and multiple imputation (MI)—and their adaptation to geostatistical workflows.

2. Multivariate transforms

This section introduces required notation and multivariate stationary random functions. Next, the section describes the use of the LMC to represent the spatial correlation structure of coregionalized variables. Fitting an LMC with sparse multivariate data is difficult in practice so decorrelation transforms have become a common alternative; these remove correlation between variables so the model can be generated independently for each factor and recombined to model the target variables. An overview of the most popular decorrelation transformations is provided. The section concludes with a discussion of practical considerations and implementation details of the multivariate data transformations.

2.1 Multivariate stationary random functions

Stationarity is an important concept in modeling of continuous spatial variables. Consider that $\left\{\mathbf{Z}\left(\mathbf{u}\right)\mathbf{:}\mathbf{u}\in D\right\}$ is a second-order stationary random function defined over a region D and $\left\{{\mathrm{z}}_{\mathrm{i}}\left({\mathbf{u}}_k\right),k=1,\cdots ,n;i=1,\cdots ,m\right\}$ represents n observations of m correlated variables. A random function is second-order stationary when the mean vector of $\mathbf{Z}\left(\mathbf{u}\right)$ is calculated by:

$$E\left[\mathbf{Z}\left(\mathbf{u}\right)\right]=\mathbf{m},\forall \mathbf{u}\in D$$E1

and the variance vector ${\mathit{\sigma }}^{\mathbf{2}}$ is defined as:

$$E\left[{\left(\mathbf{Z}\left(\mathbf{u}\right)-\mathbf{m}\right)}^2\right]={\mathbf{\sigma }}^{\mathbf{2}},\forall \mathbf{u}\in D$$E2

The variance-covariance matrix at a lag $\mathbf{h}=0$ is $\mathbf{C}\left(0\right)=\mathit{cov}\left(\mathbf{Z}\left(\mathbf{u}\right)\right)$ and the covariance is:

$$E\left[{\left(\mathbf{Z}\left(\mathbf{u}\right)-\mathbf{m}\right)}^T\left(\mathbf{Z}\left(\mathbf{u}+\mathbf{h}\right)-\mathbf{m}\right)\right]=\mathbf{C}\left(\mathbf{h}\right),\forall \mathbf{u},\mathbf{u}+\mathbf{h}\in D$$E3

The variogram matrix of direct and cross-variograms is given by:

$$\mathbf{\Gamma }\left(\mathbf{h}\right)=\left(1/2\right)E\left[{\left(\mathbf{Z}\left(\mathbf{u}\right)-\mathbf{Z}\left(\mathbf{u}+\mathbf{h}\right)\right)}^T\left(\mathbf{Z}\left(\mathbf{u}\right)-\mathbf{Z}(\mathbf{u}+\mathbf{h}\right)\right],\forall \mathbf{u},\mathbf{u}+\mathbf{h}\in D$$E4

where $T$ represents transposition, $\mathbf{Z}\left(\mathbf{u}\right)=\left[Z_1\left(\mathbf{u}\right),\dots ,Z_m\left(\mathbf{u}\right)\right]$, and $\mathbf{m}=\left[m_1,\dots ,m_m\right]$. If a random function is second-order stationary and also intrinsic, the variogram $\mathbf{\Gamma }\left(\mathbf{h}\right)$ and its covariance $\mathbf{C}\left(\mathbf{h}\right)$ are equivalent functions; that is, $\mathbf{\Gamma }\left(\mathbf{h}\right)$=$\mathbf{C}\left(\mathbf{0}\right)\mathbf{-}\mathbf{C}\left(\mathbf{h}\right)$.

2.2 The linear model of coregionalization (LMC)

LMC assumes that $Z_i\left(\mathbf{u}\right)$ is a linear combination of pairwise and spatially uncorrelated factors $Y^k\left(\mathbf{u}\right)$, $k=0,\dots ,K$ with $E[Y^k\left(\mathbf{u}\right)$] = 0 and $\mathit{var}(Y^k\left(\mathbf{u}\right)$) = 1 [18, 19, 20]; that is,

$$Z_i\left(\mathbf{u}\right)={\mathrm{m}}_{\mathrm{i}}+\sum _{k=0}^K{\alpha }_i^k{Y}_i^k\left(\mathbf{u}\right),\forall i,\mathbf{u}$$E5

An LMC fits direct and cross-variograms jointly with a positive definite model. In this procedure, $Y_i^k\left(\mathbf{u}\right)$ is assumed to have the same variogram function ${\gamma }^k\left(\mathbf{h}\right)$ with $b_{\mathit{ij}}^k={\alpha }_i^k{\alpha }_j^k$, $\forall i,j$ coefficients, with the matrix form of the LMC given by:

$$\mathbf{\Gamma }\left(\mathbf{h}\right)=\sum _{k=0}^K{\mathbf{B}}^k{\gamma }^k\left(\mathbf{h}\right)$$E6

where ${\mathbf{B}}^k$ is a positive semi-definite coregionalization matrix. This matrix contains the contributions of the variogram functions ${\gamma }^k\left(\mathbf{h}\right)$. First the $K+1$ basic variogram functions ${\gamma }^k\left(\mathbf{h}\right)$ are defined with common range parameters for all direct and cross-variograms of the individual variogram functions. An automatic fitting procedure can be implemented to iteratively determine the coefficients $b_{\mathit{ij}}^k$ that fit the experimental data [21].

2.3 Principal component analysis (PCA)

PCA is a process that transforms the set of variables to a new reference where they are uncorrelated (Figure 1) [5, 18]. PCA serves many purposes in geostatistical modeling such as dimensionality reduction, identification of underlying factors, decorrelation of variables for independent modeling of multivariate systems, and verification of intrinsic correlation [22].

Figure 1.

Original coordinates/variables (left) and rotated uncorrelated coordinates/variables along the principal axis of the data through PCA (right).

Consider the correlated variables $\mathbf{Z}\left(\mathbf{u}\right)=\left[Z_1\left(\mathbf{u}\right),\dots ,Z_m\left(\mathbf{u}\right)\right]$ are transformed into pairwise ($\mathbf{h}=0$) uncorrelated factors $\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)=\left[{\stackrel{\sim }{Y}}_1\left(\mathbf{u}\right),\dots ,{\stackrel{\sim }{Y}}_m\left(\mathbf{u}\right)\right]$. Figure 1 (right) shows two orthogonal uncorrelated PCA factors (principal components). ${\stackrel{\sim }{Y}}_1$ accounts for the majority of variance within the data, while ${\stackrel{\sim }{Y}}_2$ accounts for the remaining variance not explained by ${\stackrel{\sim }{Y}}_1$. The implementation of the PCA transform and its back-transform are given as follows [23]:

1. Transformation of correlated variables $\mathbf{Z}\left(\mathbf{u}\right)=\left[Z_1\left(\mathbf{u}\right),\dots ,Z_m\left(\mathbf{u}\right)\right]$ to Gaussian variables $\mathbf{Y}\left(\mathbf{u}\right)\sim N\left(0,1\right)$. $\mathbf{Y}\left(\mathbf{u}\right)={\mathbf{G}}^{\mathbf{-}\mathbf{1}}\left[\mathbf{F}\left(\mathbf{Z}\left(\mathbf{u}\right)\right)\right]$. This step is optional but having a mean of zero and a variance of one improves the interpretability of the results.

2. Calculation of the covariance matrix of $\mathbf{Y}\left(\mathbf{u}\right)$; that is, ${\sum }_{\mathrm{Y}}=\left(\frac{1}{n}\right)\left[\mathbf{Y}\left(\mathbf{u}\right)\mathbf{Y}{\left(\mathbf{u}\right)}^T\right]\mathit{,}$

   where Σ_Y symbolizes the variance-covariance matrix.

3. Spectral decomposition of ${\sum }_{\mathbf{Y}}$ yields the eigenvector matrix V$\left({\mathbf{V}}^{\mathbf{T}}\mathbf{V}\mathbf{=}\mathbf{I}\right)$ and the diagonal eigenvalue matrix D, ${\sum }_{\mathbf{Y}}=\mathbf{VD}{\mathbf{V}}^T$. While each row of the orthogonal matrix $\mathbf{V}$ has vectors that explain the primary directions of the variability, the magnitude of the variability for each row of $\mathbf{V}$ is described by the diagonal matrix $\mathbf{D}$.

4. Calculation of the standardized components $\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$. This is the linear matrix multiplication of $\mathbf{Y}\left(\mathbf{u}\right)$ with $\mathbf{B}\mathbf{,}\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)\mathbf{=}\mathbf{BY}\left(\mathbf{u}\right)$, where $\mathbf{B}\mathbf{=}{\mathbf{D}}^{\frac{\mathbf{-}\mathbf{1}}{\mathbf{2}}}.$

5. Back transformation of $\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$ to the Gaussian variables; that is, $\mathbf{Y}\left(\mathbf{u}\right)\mathbf{=}{\mathbf{V}}^{\mathbf{T}}\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$, if used in the modeling workflow.

6. Back transformation of $\mathbf{Y}\left(\mathbf{u}\right)$ to the original variables $\mathbf{Z}\left(\mathbf{u}\right)$; that is, $\mathbf{Z}\left(\mathbf{u}\right)\mathbf{=}{\mathbf{F}}^{\mathbf{-}\mathbf{1}}\left[\mathbf{G}\left(\mathbf{Y}\left(\mathbf{u}\right)\right)\right]$.

The PCA transform and back transform are schematically illustrated in Figure 2 with the Jura dataset [24, 25]. Note that only two variables (nickel (Ni) and zinc (Zn)) are included for clarity.

Figure 2.

Schematic representation of the PCA transform and back-transform steps for Ni and Zn variables of the Jura dataset.

The correlated Ni and Zn variables of Jura data in Figure 2 are first transformed into normal scores (NS_Ni and NS_Zn). Next, the components are generated (${\stackrel{\sim }{Y}}_1$ and ${\stackrel{\sim }{Y}}_2$) for the data. Finally, the components are back-transformed to normal scores and then a normal score back-transform is performed to regenerate the original values. The principal components (PCs) are obtained in a ranked order, where the first component is given by the first eigenvector of the spectral decomposition of ${\sum }_{\mathbf{Y}}$ and is associated with the first eigenvalue. The eigenvalues are linked to the proportion of variance associated with each PC. Generally, the majority of data variance is accounted for within the first few PCs, with the proportion calculated by the following equation:

$$d_{i,i}/\mathit{tr}\left(\mathbf{D}\right)\bullet 100$$E7

$d_{i,i}$ describes each diagonal entry in $\mathbf{D}$ that corresponds to the variance of a PC. The eigenvalues of each component of the Jura data, for example, in Figure 2 are 1.66 (accounting for 85% of the total variance) and 0.34 (accounting for 15% of the total variance), respectively. By obtaining the proportion of data variance represented, the relative variance contribution of each component can be identified. Knowing each component’s variance contribution provides insight into whether a dimensionality reduction is necessary to reduce system complexity. Dimensionality reduction retains components with high system-variance contributions. This simplifies the overall system while avoiding significant information loss. However, it is important to highlight that PCA is sensitive to outliers. In the presence of outliers, the ratio of variance attributed to each component may be distorted. Thus, it is recommended that the data be standardized before transformation. Moreover, if fewer than m components are used in geostatistical analysis, the back transformation will not perfectly reproduce the input data, which is typically required in spatial modeling contexts.

PCA is widely used due to its straightforward application. Nonetheless, there are drawbacks when the goal is to remove spatial correlation from the data because decorrelation can only be ensured for null lag separation vectors ($\mathbf{h}=0)$. PCA can only achieve spatial decorrelation in intrinsic correlation models where the correlation between variables is independent of the lag$\mathbf{h}$ [26]. PCA analysis is rather useful for detecting intrinsic correlation among variables. It can be directly applied in multivariate modeling procedures when the correlation between two variables does not depend on the spatial scale. In such cases, the direct and cross-covariance functions are chosen proportional to the same basic spatial correlation function. On the other hand, when the variables do not conform to an intrinsic correlation model, the use of alternative transformations tailored for spatial decorrelation should be considered.

PCA also assumes that the $\mathbf{h}=0$ multivariate distribution is completely parameterized by ${\sum }_{\mathbf{Y}}$ [13]. When the data has complex features, such as nonlinear relationships, mineralogical constraints, and heteroscedasticity, this assumption causes uncorrelated but not independent principal components at $\mathbf{h}=0$. The back transformation of those components is unlikely to reproduce the original multivariate dependencies.

2.4 Minimum/maximum autocorrelation factors

The MAF transform was first proposed by Switzer [8] for use in image processing and remote sensing fields. It was later used by ref. [7] in geostatistical modeling. MAF applies PCA transforms twice, to the variance-covariance matrix ${\sum }_{\mathbf{Y}}$ and the experimental variogram matrix of the standardized components ${\widehat{\mathbf{\Gamma }}}_{\stackrel{\sim }{\mathbf{Y}}}\left(\mathbf{h}\right)$ built on one chosen lag $\left(\mathbf{h}\ne 0\right).$ The MAF transformation results in random variables that are spatially uncorrelated at $\mathbf{h}=0$ and have no cross-correlations at the chosen lag ${\mathbf{h}}_j>0,j=1,\dots ,J$ [27, 28]. The implementation steps of the MAF transform and its back-transform are given as follows:

1. Calculate the standardized components $\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$ at $\mathbf{h}=0$. Refer to PCA steps I-IV given in Section 2.3.

2. Select the lag distance $\left\{{\mathbf{h}}_j,j=1,\dots ,J\right\}$ and compute the experimental variogram matrix $\widehat{\mathbf{\Gamma }}\left({\mathbf{h}}_j\right)$.

3. Diagonalize $\widehat{\mathbf{\Gamma }}\left({\mathbf{h}}_j\right)$; that is $\widehat{\mathbf{\Gamma }}\left({\mathbf{h}}_j\right)$=${\mathbf{V}}_{\mathrm{Y}}\mathbf{D}{\mathbf{V}}_{\mathrm{Y}}^{\mathrm{T}}$.

   Calculate the spatially uncorrelated MAF factors ${\stackrel{\sim }{\mathbf{Y}}}^{\mathbf{\ast }}\left(\mathbf{u}\right)$; that is, ${\stackrel{\sim }{\mathbf{Y}}}^{\mathbf{\ast }}\left(\mathbf{u}\right)\mathbf{=}{\mathbf{V}}_{\mathrm{Y}}\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$. This step removes structure at both $\mathbf{h}=0$ and ${\mathbf{h}}_j>$0

4. Back transform ${\stackrel{\sim }{\mathbf{Y}}}^{\mathbf{\ast }}\left(\mathbf{u}\right)$ to the PCs; that is,$\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)\mathbf{=}{\mathbf{V}}_{\mathbf{Y}}^{\mathbf{-}\mathbf{1}}{\stackrel{\sim }{\mathbf{Y}}}^{\mathbf{\ast }}\left(\mathbf{u}\right)$.

5. Back transform $\stackrel{\sim }{\mathbf{Y}}\left(\mathbf{u}\right)$ to the Gaussian $\mathbf{Y}\left(\mathbf{u}\right)$ and the original values$\mathbf{Z}\left(\mathbf{u}\right)$. Refer to steps V and VI in Section 2.3.

This two-step transformation generates factors ${\stackrel{\sim }{\mathbf{Y}}}^{\mathbf{\ast }}\left(\mathbf{u}\right)$ that are uncorrelated at both $\mathbf{h}=0$ and the selected lag distance ${\mathbf{h}}_j>0$ (steps III, IV). Figure 3 shows the scatter plots of the original chromium (Cr), Ni, and Zn variables of the Jura dataset and the scatter plots of the MAF factors.

Figure 3.

Scatter plots of the correlated data (left) and MAF data (right) for Cr, Ni, and Zn variables of the Jura dataset.

It is clear from Figure 3 that the MAF transform produces uncorrelated factors MAF1, MAF2, and MAF3. The important aspect of the MAF transformation is found in the cross-variograms. Figure 4 shows the cross-variograms of the original and transformed variables. Although MAF is efficient and facilitates a simplified modeling process as the factors are uncorrelated, there are limitations. For example, spatial decorrelation only occurs for a specific lag distance $\mathbf{h}$ that is user-defined. The choice of $\mathbf{h}$ is a modeling decision, but it is commonly chosen according to the lag in which the spatial correlation is observed for components (i.e., PCAs) obtained in a previous exploratory data analysis step.

Figure 4.

The original and transformed variables cross-variograms for Cr, Ni, and Zn variables of the Jura dataset.

The cross-correlations for MAF in Figure 4 are zero at $h=0.6\mathit{km}$. Dimitrakopoulos and Fonseca [29] presented a model-based derivation of MAF, where the direct and cross-spatial continuity of the attributes are used to transform the original variables and obtain uncorrelated factors at any lag distance. Vargas-Guzmán [30] discuss that MAF factors are generally based on omnidirectional covariance matrices, but each MAF factor can be modeled according to its anisotropy; however, the solution may not be unique if direction is considered [30]. Bailey and Krzanowski [31] proposed a generalization of MAF that produces factors with nearly zero cross-correlations at all given lags using a rotation matrix. Nonetheless, like PCA, MAF is a linear transformation that does not account for complex multivariate features, which implies that such relations may not be reproduced in the final model. This has motivated the use of multivariate transforms such as SCT and PPMT.

2.5 Stepwise conditional transform

Many geostatistical workflows rely on multivariate Gaussian distributions; however, geological variables are rarely Gaussian. For this reason, variables are often transformed prior to geostatistical modeling [12]. Generally, normal scores transform [31, 32] is applied, but this approach only ensures the transformed variables have univariate Gaussian distributions. In the presence of multivariate correlated variables, the modeling process becomes increasingly difficult. Independence between variables cannot be guaranteed using linear decorrelation transforms (PCA and MAF) in datasets with complex multivariate features [12]. SCT [9] and PPMT [12] are alternative transformations that aim to generate components that are independent, rather than just uncorrelated.

SCT was first introduced by Rosenblatt [10] and applied to geostatistical workflows by Leuangthong and Deutsch [9]. This transform handles complex features while decorrelating the variables at $\mathbf{h}=0$ to create uncorrelated independent multiGaussian distributions. SCT is performed in a chained manner, where the $\mathit{Mth}$ variable is conditionally transformed based on the first $\left(M-1\right)$ variables [9]. The steps in stepwise conditional transformation are given as follows:

1. Transformation of the first variable $Z_1\left(\mathbf{u}\right)$ to a Gaussian variable $Y_1\left(\mathbf{u}\right)\sim N\left(0,1\right)$; that is, $Y_1\left(\mathbf{u}\right)={\mathbf{G}}^{\mathbf{-}\mathbf{1}}\left[{\mathbf{F}}_{\mathbf{1}}\left({\mathbf{Z}}_{\mathbf{1}}\left(\mathbf{u}\right)\right)\right]$.

2. Conditional transformation of subsequent variables $Z_2\left(\mathbf{u}\right)\cdots Z_M\left(\mathbf{u}\right)$ to the previous variables $\left(M-1\right)$; that is, $Y_M\left(\mathbf{u}\right)={\mathbf{G}}^{\mathbf{-}\mathbf{1}}\left[{\mathbf{F}}_{\mathbf{M}\mid \mathbf{1}\mathbf{,}\dots \mathbf{,}\mathbf{M}\mathbf{-}\mathbf{1}}\left({\mathbf{Z}}_{\mathbf{M}}\left(\mathbf{u}\right)|{\mathbf{Y}}_{\mathbf{1}\mathbf{,}\mathbf{\cdots }\mathbf{,}\mathbf{M}\mathbf{-}\mathbf{1}}\left(\mathbf{u}\right)\right)\right]$.

3. Back transformation of the $\mathit{Mth}$ Gaussian variable $Y_M\left(\mathbf{u}\right)$ to the original variables $Z_M\left(\mathbf{u}\right)$ that is, $Z_M\left(\mathbf{u}\right)={\mathbf{F}}_{\mathrm{M}\mid 1,\dots ,\mathrm{M}-1}^{\mathbf{-}\mathbf{1}}\left[\mathbf{G}\left(Y_M\left(\mathbf{u}\right)|{\mathrm{Y}}_{1,\cdots ,\mathrm{M}-1}\left(\mathbf{u}\right)\right)\right]$.

SCT facilitates the simulation of multiple variables with complex relationships by enabling the independent Gaussian simulation of $Y_M\left(\mathbf{u}\right)$ . The SCT back transform provides the reproduction of the original complex features. This procedure is illustrated in Figure 5 for the correlated variables (Ni and Zn) from the Jura dataset.

Figure 5.

SCT and its back-transform for Ni and Zn variables of the Jura dataset.

In Figure 5, a normal score transformation is applied to Ni, while Zn is conditionally transformed based on Ni. For SCT, the order of the variables in the transformation affects the final components obtained; typically, the most sampled or the most important variables are transformed first. The first variable does not require any other variables to define the transform, while the subsequent variables are conditional to all previously transformed variables and their variograms are the integration of the spatial structure of previously transformed variables and their cross-correlation [9].

Another challenge with SCT is that many samples are required to obtain reliable conditional distributions. Although Leuangthong and Deutsch [9] proposed that the gaps in the multivariate distribution can be filled through kernel density estimation [33], SCT may not be feasible for more than two to four variables in geological datasets because of typical sample quantities [13]. SCT is also prone to artifacts at the binning boundaries. Thus, fitting a probabilistic Gaussian mixture model (GMM) can be useful to calculate the conditional distributions required for SCT and avoid binning [34].

SCT has been successfully applied in geostatistical workflows. For example, Hosseini and Asghari [35] applied SCT to the simulation of three variables with complex multivariate relationships with stoichiometric constraints. Recently, Figueiredo et al. [36] proposed an algorithm that performs direct simulation using the SCT. The algorithm introduces a uniformly distributed auxiliary random variable that establishes conditional distributions non-parametrically. These distributions are then used to form the spatial correlation model.

2.6 Projection pursuit multivariate transform (PPMT)

The projection pursuit algorithm was initially proposed by Friedman and Tukey [37] for multivariate data analysis. The core idea of this algorithm is to find directions onto which multidimensional data can be projected to reveal important data features. Projection pursuit has been applied to density estimation, regression, and classification [38]. The algorithm also has a factor that transforms the original data toward a multivariate Gaussian distribution [37, 39]; PPMT was proposed to make use of this tendency within a geostatistical framework [12]. PPMT transforms complex high-dimensional geologic data to a pairwise ($\mathbf{h}=0)$ uncorrelated multi-Gaussian distribution. It is a non-linear extension of PCA that consists of two steps: (1) preprocessing to remove the linear dependency between the variables and (2) using projection pursuit to iteratively remove complex multivariate features and create independent variables. The implementation steps of the PPMT forward transform and its back-transform are given as follows:

1. Preprocess the data: Follow steps I–III given in Section 2.3, considering the Gaussian transform.

2. Transform the data $\mathbf{Y}\left(\mathbf{u}\right)$ to be uncorrelated with unit variance by sphering; that is ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)\mathbf{=}{\mathbf{S}}^{\frac{\mathbf{-}\mathbf{1}}{\mathbf{2}}}\mathbf{Y}\left(\mathbf{u}\right)$, where ${\mathbf{S}}^{\frac{\mathbf{-}\mathbf{1}}{\mathbf{2}}}\mathbf{=}\mathbf{V}{\mathbf{D}}^{\frac{\mathbf{-}\mathbf{1}}{\mathbf{2}}}{\mathbf{V}}^{\mathbf{T}}$. The rotated data ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$ has an identity covariance matrix.

   These preprocessing steps are visualized for the correlated variables (Cr and Zn) from the Jura dataset in Figure 6.

   It is clear from Figure 6 that the sphere transform removed the bivariate relationship in the original and normal score data. Although sphering only removes the linear correlation, complex dependencies still exist between the variables. The PPMT addresses these complexities. After preprocessing, the remaining PPMT steps are:

3. Consider an $m$ dimension unit length vector $\mathit{\theta }$ and the projection of ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$ upon it, $\mathbf{p}\mathbf{=}\mathit{\theta }{\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$.

4. Determine the optimum $\mathit{\theta }$ that maximizes the projection index $\mathit{I}\left(\mathit{\theta }\right)$ quantifying its univariate non-Gaussian nature. If $\mathbf{p}$ is perfectly Gaussian for any $\mathit{\theta }$, then $\mathit{I}\left(\mathit{\theta }\right)$ is equal to zero. In this step, ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$ is transformed into $\widehat{\mathbf{Y}}\left(\mathbf{u}\right)$ so that the projection $\mathbf{p}$ along $\mathit{\theta }\mathbf{,}\widehat{\mathbf{p}}\mathbf{=}\mathit{\theta }\widehat{\mathbf{Y}}\left(\mathbf{u}\right)$ is Gaussianized [40].

5. Iterate Step 4 until $\mathit{I}\left(\mathit{\theta }\right)$ reaches convergence. This yields an uncorrelated and multivariate Gaussian $\widehat{\mathbf{Y}}\left(\mathbf{u}\right)$. Studies discussed in Refs. [12, 40] highlight that defining the target value for $\mathit{I}\left(\mathit{\theta }\right)$ is complex. For the details on the stopping criteria, the interested reader is referred to Ref. [37].

6. Reverse of the forward transform order to back-transform $\widehat{\mathbf{Y}}\left(\mathbf{u}\right)$ to ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$, for each projection-pursuit iteration [40].

7. Back-transform ${\mathbf{Y}}^{\prime }\left(\mathbf{u}\right)$ to the Gaussian variables $\mathbf{Y}\left(\mathbf{u}\right)$; that is, $\mathbf{Y}\left(\mathbf{u}\right)\mathbf{=}{\mathbf{Y}}^{\prime }\left(\mathbf{u}\right){\mathbf{S}}^{\frac{\mathbf{1}}{\mathbf{2}}}$.

8. Back-transform $\mathbf{Y}\left(\mathbf{u}\right)$ to the original variables $\mathbf{Z}\left(\mathbf{u}\right)$; that is, $Z\left(\mathbf{u}\right)\mathbf{=}{\mathbf{F}}^{\mathbf{-}\mathbf{1}}\left[\mathbf{G}\left(\mathbf{Y}\left(\mathbf{u}\right)\right)\right]$.

Figure 6.

The original, normal score, and sphered data scatter plots.

The PPMT transform and the convergence of the projection index are illustrated with the Jura data in Figure 7.

Figure 7.

Scatter plots of the PPMT transform (left) and progression of the projection index (right) over 65 iterations.

Figure 7 (right) shows the convergence of the projection index $\mathit{I}\left(\mathit{\theta }\right)$ through the iteration process. The non-Gaussian nature of the projection decreases (i.e., it becomes more Gaussian-like) after 60 iterations. Increase in the projection index across iterations indicate that a local optimum is found in the previous iteration.

PPMT can be applied to any number of variables. However, similar to PCA and SCT, it obtains decorrelation of the variables at a null $\left(\mathbf{h}=0\right)$ lag distance. Thus, the spatial dependency must be inspected after the transform is performed. In the case of residual spatial correlation between the transformed variables, MAF may be applied to remove remnant correlations at the lag distance of interest. Kloeckner et al. [41] proposed that the modeling effort in multivariate settings can be greatly reduced by using the PPMT where each transformed variable (Step V. above) can be modeled independently. Kloeckner et al. [41] showed that the transformation process did not cause significant distortions in the data distributions and that key statistical parameters, such as histograms and spatial correlations, were reproduced in the final model.

2.7 Summary

In geostatistical modeling, the use of multivariate transformations is an effective alternative to existing intensive multivariate modeling workflows like LMC. The multivariate transformations apply a decorrelation or independence transformation, allowing each resultant factor to be modeled independently. The first and second-order statistics and original bivariate relationships among the data are restored through the back transform. PCA (see Section 2.3) decomposes the correlated multivariate variables into orthogonal factors that are pairwise uncorrelated at a lag distance zero $\left(\mathbf{h}=0\right)$. PCA assumes that the factors are also uncorrelated at all possible lags $\left\{{\mathbf{h}}_j,j=1,\dots ,J\right\}$ - though only valid in the case of intrinsic correlation models. MAF (see Section 2.4) is a linear extension of PCA and transforms variables into orthogonal factors; these are uncorrelated at $\mathbf{h}=0$ and have no cross-correlations at a particular lag distance ${\mathbf{h}}_{\mathrm{j}}>0$. MAF factors are uncorrelated at all lag distances assuming that the multivariate data can be modeled using a two-structure LMC. As the PCA and MAF are linear transforms, they do not result in independent variables in cases with multivariate complexities such as non-linearity, heteroscedasticity, and constraints. SCT (Section 2.5) and PPMT (Section 2.6) are the most common multivariate Gaussian transforms and they remove non-Gaussian features, resulting in transformed multi-Gaussian data that are independent at the lag distance zero $\left(\mathbf{h}=0\right)$. While the implementation of SCT is often challenging for more than two to four variables, PPMT can manage datasets containing a high number of variables. Moreover, the combination of PPMT and MAF can be used as a chain transform when correlation is observed at a lag >0.

3. Data imputation

In a mineral deposit, it is common to have missing variables at sampling locations; sampling cost, changing field conditions, sampling quality, evolving drilling technology, and changing project goals mean that the variables assayed change over time and space. Unequally sampled (heterotopic) data presents a problem in geostatistical analysis when the selected multivariate transformation (Section 2) requires equally sampled (homotopic) data at all data locations. Figure 8 demonstrates this missing data concern for a simple example with heterotopic data.

Figure 8.

A dataset composed of three heterotopically sampled variables at nine locations. There are only 2 locations (4 and 7) where all three variables are collocated. Samples in different colors on the map are shown in the figure as slightly offset for clarity, but they are collocated.

The sample data in Figure 8 contains 18 data entries at 9 locations and three variables of interest. The blank data entries indicate missing variables in a particular sample. In this example, multivariate transformation methods cannot be applied directly to any record where one or more variables are missing. One option is to exclude incomplete observations; however, this creates a loss of information that may introduce additional uncertainty and bias into the resultant model [14]. To maintain the integrity of the original data, imputation of the missing variables through estimation or simulation can be considered as an alternative approach to satisfy a homotopic data requirement in the modeling workflow.

In the following section, missing data mechanisms that cause data to be missing are introduced. This is followed by an overview of the conventional missing data imputation methods and a discussion of their suitability in various missing data mechanisms. Imputation methods that consider the spatial and multivariate structure of geological data are reviewed. The section concludes with a discussion regarding the application of imputation methods within geostatistical workflows.

3.1 Missing data mechanisms

Some of the earliest considerations of missing data mechanisms were provided by Rubin [42], and the reader is referred there for more details. The more general ideas are discussed here.

Consider that K spatial variables $Z_1\dots Z_K$ have been sampled in differing subsets at $n$ observation locations, resulting in the data matrix $\mathbf{Z}$: $z_{i,k},i=1,\dots ,n,k=1,\dots ,K$. According to Rubin [42], a binary matrix $\mathbf{R}$ can be defined to specify the sampled (0) and missing values (1) of $\mathbf{Z}$ associated with a probability distribution $f\left(\mathbf{R}|\mathbf{Z},\varphi \right),$where $f$ is the probability distribution function and $\varphi$is a vector of parameters that explain the missingness.

Missingness is grouped into three categories, arising from different mechanisms: (1) missing completely at random (MCAR), (2) missing at random (MAR), and (3) missing not at random (MNAR) [43]. The MCAR missing mechanism occurs when the probability of being missing is independent of the values $\mathbf{Z}$ according to $f\left(\mathbf{R}|\mathbf{Z},\varphi \right)=f\left(\mathbf{R}|\varphi \right)\forall \mathbf{Z}$,$\varphi$. Consider that ${\mathbf{Z}}_{\mathrm{O}}$ and ${\mathbf{Z}}_{\mathrm{M}}$ represent the observed and missing values, respectively. MCAR implies that ${\mathbf{Z}}_{\mathrm{M}}$ and ${\mathbf{Z}}_{\mathrm{O}}$ share the same overall distributions. This mechanism does not create a bias in further statistical analysis and geostatistical modeling.

MAR occurs when the probability of missingness is only related to ${\mathbf{Z}}_{\mathrm{O}}$ according to $f\left(\mathbf{R}|\mathbf{Z},\varphi \right)=f\left(\mathbf{R}|{\mathbf{Z}}_{\mathrm{O}},\varphi \right)\forall {\mathbf{Z}}_{\mathrm{M}},\varphi$. MAR implies that missingness is random for ${\mathbf{Z}}_{\mathrm{M}}$. MAR is also referred to as ignorable missingness. Available imputation methods may generate unbiased results without special handling in the case of MAR.

If missing mechanisms are neither MCAR nor MAR, then the probability of being missing is not random and is related to ${\mathbf{Z}}_{\mathrm{O}}$ and ${\mathbf{Z}}_{\mathrm{M}}$ according to $f\left(\mathbf{R}|\mathbf{Z},\varphi \right)=f\left(\mathbf{R}|{\mathbf{Z}}_{\mathrm{O}},{\mathbf{Z}}_{\mathrm{M}},\varphi \right)\forall \varphi$ [44]. The MNAR mechanism is the most complex; the missing data hold important information about the data population that cannot be ignored and introduces bias in subsequent modeling [45]. Figure 9 illustrates the three different missing data mechanisms (MCAR, MAR, and MNAR) with the Walker Lake dataset [46]. The Walker Lake dataset used herein consists of two correlated variables (U and V, which were anonymously named, representing concentration in parts per million).

Figure 9.

Walker Lake dataset example depicting the three missing data mechanisms (MCAR, MAR, and MNAR) for the variable U.

In Figure 9, the complete dataset for the variable U was produced by regularly sampling the dataset at 3068 locations. The MCAR mechanism was created by setting randomly drawn U values to be missing. V values below 135.49 and U values below 20.49 were assigned as missing for the MAR and MNAR mechanisms, respectively. The range of missing values was set to 35% for each missing mechanism. Figure 9 also shows how missing values in different mechanisms affect the U variable distribution. U variable statistics do not change significantly in the MCAR mechanism. In comparison, there are variations in the statistics of U for both MAR and MNAR. Figure 10 illustrates the effect of different missing mechanisms on the variogram maps.

Figure 10.

Walker Lake dataset variogram maps depicting changes in anisotropy direction using various missing data mechanisms for the variable U.

It is clear from Figure 10 that while the orientation of MCAR appears close to the observed variogram map of the complete data, the orientation of MAR and MNAR differ because of the mechanism of missing data.

3.2 Imputation methods

Imputation can be defined as the replacement/modeling of missing values with a suitable method. Many imputation methods exist in literature [39]; the most common include single imputation (SI), maximum likelihood estimation (MLE), and multiple imputation (MI). An overview and discussion of these three methods are given below.

3.2.1 Single imputation (SI)

These methods generate a single value to replace missing data based on a certain mathematical rule. Single imputation methods include the following:

1. Arithmetic mean imputation: This method fills in missing values with the stationary mean of that variable [47]. This method alters the variance of the distribution but preserves the mean. This method may not be effective in mineral deposit applications where cut-off grades are important.

2. Regression or conditional mean imputation: This method fills in missing values based on a regression fit to correlated secondary variable(s) available at all missing locations [48]. Cokriging could be used. This method may produce biased parameter estimates depending on the missing mechanism and relationship to the secondary variable.

3. Stochastic regression imputation: This method is an improvement over regression imputation. It applies a stochastic method to add realistic variability to the regression model under the assumption of normally distributed errors [49]. This method removes the biases in standard regression imputation methods.

Figure 11 shows the application of single imputation methods on a simplified example with a MAR mechanism. Figure 11(a) displays the complete dataset, with gray points representing the observed data and red points indicating the missing data that were initially removed from the overall dataset to demonstrate the single SI methods.

Figure 11.

Single imputation methods: (a) complete dataset, red points are removed in subsequent plots for demonstration, (b) mean imputation of the red points, (c) regression imputation of the red points, (d) stochastic regression imputation of the red points.

SI methods, apart from stochastic regression imputation, may introduce bias for MAR and MNAR mechanisms. Stochastic regression imputation typically results in less bias. SI imputation does not consider uncertainty in imputed values. This is the main motivation for other methods.

3.2.2 Maximum likelihood estimation (MLE)

MLE uses a log-likelihood function to identify distribution model parameters $\left(\theta \in R\right)$, such as mean and variance, with the greatest likelihood of producing the observed data. A general method to find the unknown population $\theta$ is the expectation-maximization (EM) algorithm, this is useful when the joint distribution of ${\mathbf{Z}}_{\mathrm{O}}$ and ${\mathbf{Z}}_{\mathrm{M}}$ is explicit [50]. The basic idea of EM is to identify the estimate of $\theta$ that maximizes the observed data log-likelihood $l\left(\theta ,{\mathbf{Z}}_{\mathrm{O}}\right)$, that is, the probability density function of the observations $f\left({\mathbf{Z}}_{\mathrm{O}};\theta \right)$:

$$l\left(\theta ,{\mathbf{Z}}_{\mathrm{O}}\right)=\mathit{\text{logf}}\left({\mathbf{Z}}_{\mathrm{O}};\theta \right)=\mathit{log}\int f\left({\mathbf{Z}}_{\mathrm{O}},{\mathbf{Z}}_{\mathrm{M}};\theta \right)d{\mathbf{Z}}_{\mathrm{M}}$$E8

In general cases, the quantity described in Eq. 8 cannot be calculated explicitly. Therefore, EM defines the MLE by iteratively maximizing the expected complete-data log-likelihood $l\left(\mathbf{Z},\theta \right)$ as follows:

$$l\left(\mathbf{Z},\theta \right)=\mathit{\text{logf}}\left({\mathbf{Z}}_{\mathrm{O}},{\mathbf{Z}}_{\mathrm{M}};\theta \right)$$E9

Consider that ${\theta }^{\left(0\right)}$ represents an arbitrary initial value and ${\theta }^{\left(t\right)}$ is the estimate of $\theta$ at $t^{\mathit{th}}$ iteration. The EM algorithm implements two steps in its next iteration:

1. The expectation ($E$) step: The expectation of $l\left(\mathbf{Z},\theta \right)$ regarding the conditional distribution of missing covariate parameterized by ${\theta }^{\left(t\right)}:$

   $$Q\left(\theta ,{\theta }^{\left(t\right)}\right)=E\left[l\left(\mathbf{Z},\theta \right)|{\mathbf{Z}}_{\mathrm{O}};{\theta }^{\left(t\right)}\right]=\int l\left(\mathbf{Z},\theta \right)f\left({\mathbf{Z}}_{\mathrm{M}}|{\mathbf{Z}}_{\mathbf{O}};{\theta }^{\left(t\right)}\right)d{\mathbf{Z}}_{\mathrm{M}}$$E10

2. The maximization ($M$) step: ${\theta }^{\left(t+1\right)}$ is defined by maximizing the function$Q$:

$$l\left(\mathbf{Z},\theta \right)=\mathit{\text{logf}}\left({\mathbf{Z}}_{\mathrm{O}},{\mathbf{Z}}_{\mathrm{M}};\theta \right)$$E11

The iteration is terminated when successive estimates of $\theta$ are nearly identical. The ${\theta }^{\left(t+1\right)}$ that maximizes $Q\left(\theta ,{\theta }^{\left(t\right)}\right)$ is guaranteed to give an $l\left(\theta ,{\mathbf{Z}}_{\mathrm{O}}\right)$ that is greater than that provided by$,{\theta }^{\left(t\right)}$ [50].

3.2.3 Multiple imputation (MI)

In MI, multiple data sets are generated with different values imputed for each missing value at each location. Each complete dataset is then analyzed in the same fashion and could be said to represent ‘realizations of the data’. Uncertainty in the imputed values is obtained by combining the results (i.e., an ensemble) to be used in subsequent modeling or to form the overall estimates and associated standard errors.

In general, the use of MLE or MI is strongly preferred for MAR mechanisms due to their ability to accurately reproduce data distributions and provide appropriate estimates of the uncertainty for imputed data [51]. Moreover, these methods are less likely to introduce bias. Imputation for the MNAR mechanism, however, requires a specialized approach; none of the reviewed imputation methods are directly appropriate. MNAR mechanism-specific methods include selection models [52] and pattern mixture models [53, 54]. However, these methods make strong assumptions regarding the data and can introduce severe bias. Therefore, Rubin [42] and Collins et al. [55] suggested the inclusion of auxiliary variables to help explain the missingness.

For geostatistical analysis, MLE and MI are asymptotically equivalent [43] but MI is preferred because it produces data realizations that can be easily adapted to geostatistical simulation workflows. A typical MI method follows five basic steps [56]: (1) a conditional distribution $f\left({\mathbf{Z}}_{\mathrm{M}}|{\mathbf{Z}}_{\mathrm{O}}\right)$ is created to impute values using the observed data and a prior model, (2) a realization of the imputed values ${\mathbf{Z}}_M^l$ is generated by sampling imputed data stochastically from the modeled conditional distribution, (3) the procedure is repeated to produce $L$ complete datasets ${\mathbf{Z}}_n=\left[{\mathbf{Z}}_{\mathrm{O}},{\mathbf{Z}}_{\mathrm{M}}^l\right]$, $l=1,\dots ,L$, (4) statistical analysis is applied to each ${\mathbf{Z}}_l$, and (5) the results are combined. The MI process for the geostatistical modeling is illustrated in Figure 12.

Figure 12.

Schematic illustration of MI in geostatistical modeling [13].

Each imputed data realization in Figure 12 requires its own geostatistical realization. The theory of MI suggests that a standard number of realizations is around three to five [57]. More recent work, however, has recommended that MI should be performed for more realizations, say 20–100 [58]. In the geostatistical framework, 100–200 realizations are common to sample from the space of subsurface uncertainty.

3.3 Imputation for geological data

Commonly available MI methods are not directly suitable for geological data. In constructing the distribution of potential outcomes for a missing value, the following should be considered: (1) spatially correlated values; (2) collocated and correlated values of other variables; and (3) known complex multivariate features [13]. In recent years, imputation methods for geological data are commonly based on BU.

BU was proposed by Ren [17] to adapt the MI concept to a geostatistical setting. Consider a location $\mathbf{u}$ where the primary variable value is to be imputed. BU uses the primary data $\left\{y_P\left({\mathbf{u}}_i\right),i=1,\dots ,n-1\right\}$ which is not available at location $\mathbf{u}$ and the secondary data $\left\{x_j\left({\mathbf{u}}_i\right),j=1,\dots ,m,i=1,\dots ,n\right\}$ which are available in the entire domain [59]. Considering $y$ and $x_j$ are multi-Gaussian, the BU procedure is described as follows:

1. Prior distribution: Estimate the primary distribution ${\overline{y}}_P\left(\mathbf{u}\right)$ and variance ${\sigma }_P^2\left(\mathbf{u}\right)$ according to:

   $${\overline{y}}_P\left(\mathbf{u}\right)=\sum _{i=1}^n{\lambda }_i{y}_P\left({\mathbf{u}}_i\right)$$E12
   $$\sum _{i=1}^n{\lambda }_{i}C\left({\mathbf{u}}_i,{\mathbf{u}}_k\right)=C\left(\mathbf{u},{\mathbf{u}}_k\right),k=1,\dots ,n$$E13
   $${\sigma }_P^2\left(\mathbf{u}\right)=1-\sum _{i=1}^n{\lambda }_{i}C\left(\mathbf{u},{\mathbf{u}}_i\right)$$E14

   where weights ${\lambda }_i$ are calculated through the data-to-data covariance $C\left({\mathbf{u}}_i,{\mathbf{u}}_k\right)$ and the data-to-unknown covariance $C\left(\mathbf{u},{\mathbf{u}}_i\right)$ (Eq. 13). This is equivalent to using the normal equations as in simple kriging. Step 1 provides the distribution of potential outcomes according to the primary variable or variable being imputed.

2. Likelihood estimation: Define the likelihood distribution by mean ${\overline{y}}_L\left(\mathbf{u}\right)$ and variance ${\sigma }_L^2\left(\mathbf{u}\right)$ according to the following equations:

   $${\overline{y}}_L\left(\mathbf{u}\right)=\sum _{j=1}^m{\lambda }_j{x}_j\left(\mathbf{u}\right)$$E15
   $$\sum _{j=1}^m{\lambda }_j{\rho }_{\mathit{jk}}={\rho }_{\mathit{kY}},k=1,\dots ,m$$E16
   $$\underset{L}{\overset{2}{\sigma }}\left(\mathbf{u}\right)=1-\sum _{j=1}^m{\lambda }_i{\rho }_{\mathit{jY}}$$E17

   where $x_j\left(\mathbf{u}\right)$ is the secondary information at the unsampled location $\mathbf{u}$, and the weights ${\lambda }_j$ are calculated using the cross-correlation between the primary and secondary data ${\rho }_{\mathit{kY}}$, and cross-correlation between the different types of secondary data ${\rho }_{\mathit{jY}}$. This is equivalent to using the normal equations for linear least squares regression. Step 2 provides the distribution of potential outcomes according to the secondary variables.

3. Updated distribution: Merge the prior and likelihood distribution to define the updated distribution according to:

   $${\overline{y}}_U\left(\mathbf{u}\right)=\frac{{\overline{y}}_L\left(\mathbf{u}\right){\sigma }_P^2\left(\mathbf{u}\right)+{\overline{y}}_P\left(\mathbf{u}\right){\sigma }_L^2\left(\mathbf{u}\right)}{{\sigma }_P^2\left(\mathbf{u}\right)-{\sigma }_P^2\left(\mathbf{u}\right){\sigma }_L^2\left(\mathbf{u}\right)+{\sigma }_L^2\left(\mathbf{u}\right)}$$E18
   $${\sigma }_U^2\left(\mathbf{u}\right)=\frac{{\sigma }_L^2\left(\mathbf{u}\right){\sigma }_P^2\left(\mathbf{u}\right)}{{\sigma }_P^2\left(\mathbf{u}\right)-{\sigma }_P^2\left(\mathbf{u}\right){\sigma }_L^2\left(\mathbf{u}\right)+{\sigma }_L^2\left(\mathbf{u}\right)}$$E19

   where ${\overline{y}}_U\left(\mathbf{u}\right)$ is the updated mean and ${\sigma }_U^2\left(\mathbf{u}\right)$ is the updated variance. Since the posterior distribution is Gaussian and completely determined by ${\overline{y}}_U\left(\mathbf{u}\right)$ and ${\sigma }_U^2\left(\mathbf{u}\right)$, the full distribution of the missing values is obtained based on the spatially correlated primary and secondary variables. The first three steps are summarized in Figure 13.

   The likelihood variance map shown in Figure 13 is constant because the correlations ${\rho }_{\mathit{kY}}$ are derived from the data at the collocated locations and assumed to be constant in the region of interest. This leads to the variance also being constant [59].

4. Stochastically simulate: The first three steps are implemented for every missing value in a dataset. Then the data realizations are stochastically simulated by randomly sampling the (empirical) Gaussian cumulative distribution functions. It should be noted that the imputed values vary according to the degree of uncertainty in its updated distribution ${\sigma }_U^2\left(\mathbf{u}\right)$, while the sampled values remain constant (i.e., uncertainty is equal to 0) for each realization [60]. BU yields promising results when the missing mechanism is MAR.

Figure 13.

The Bayesian updating (BU) process [17].

BU is a practical imputation approach for geological data because it only requires a normal score transform that can be applied to heterotopic data. Moreover, in BU, spatial and multivariate conditioning data are incorporated and sets of variables at each location are naturally integrated. In addition, BU accounts for both spatial and collocated data sources and allows for the construction of uncertainty distributions for the missing values.

3.4 Summary

Unequally sampled data is a common occurrence in geospatial datasets. The homotopic data requirement of many multivariate modeling workflows highlights the importance of missing data imputation methodologies. Missing data imputation is the preferred approach to managing heterotopic data. The three common imputation methods are (1) SI, (2) MLE, and (3) MI. SI methods are suitable for the MCAR mechanism but can introduce bias when used for the MAR mechanism. For MAR, MI and MLE are preferred as they reproduce data variability and reflect the uncertainty of imputed values without introducing additional bias. The imputation methods reviewed cannot be applied to an MNAR mechanism. Thus, it is suggested that auxiliary variables be introduced to convert an MNAR mechanism to a MAR mechanism. For geostatistical analysis, MI is preferred to MLE as it generates several realizations of the missing values, which is consistent with many geostatistical simulation workflows that utilize realizations to account for subsurface uncertainty.

4. Conclusion

This chapter presents a detailed exploration of the complexities involved in handling multivariate data within typical geostatistical modeling workflows. It emphasizes the necessity for advanced multivariate transformations and imputation methods. The effectiveness of several transformation techniques, including PCA, MAF, SCT, and PPMT, are examined. The chapter further addresses the critical issue of missing data by exploring multiple imputation methods. These range from SI techniques to more advanced methods such as MLE and MI, considering the various missing data mechanisms.

Geospatial modeling problems continue to grow in complexity as computational resources increase, new modeling methodologies are developed, and more comprehensive multivariate sampling becomes common (drill holes, well logs, core scans, geophysical surveys, assays, etc.). Using all available data to the greatest extent possible to build the best model of a variable of interest necessitates a multivariate spatial workflow be adopted; in fact, any geostatistical problem with a single variable is now considered a ‘simple’ model. State-of-the-art geomodeling requires the incorporation of all available data and this chapter provides some of the necessary tools to enable multivariate spatial modeling of earth science variables.