Multivariate Cross-Validation and Measures of Accuracy and Precision

Cross-validation and performance measures are standard components in the evaluation of a geostatistical model. These are well established in the univariate case, but measures for multivariate geostatistical modeling have not received as much attention. In the case of a single target variable, the univariate approaches remain valid, but in the fully multivariate case where a vector of variables needs to be estimated, the evaluation needs to be based on all estimates simultaneously. An extension of cross-validation and associated performance measures to the fully multivariate case is presented and discussed for the case of regionalized compositions. The method is demonstrated by validating geostatistical models for two case studies: a sample drawn from a geochemical survey data set estimated with cokriging, and an application of direct sampling multiple-point simulation.


Introduction
Cross-validation and jackknifing are established methods for validating statistical models. In a geostatistical context, the model is based either on geostatistical estimation via kriging or on spatial simulation. The standard outputs include scatterplots and correlation coefficients of estimated against true values, and (possibly standardized) estimation errors against estimates, accompanied by error statistics (Webster and Oliver 2007). These are widely used to validate and optimize model parameters or to determine the most suitable model from a set of competing models. In addition, scatterplots of coverage probabilities versus theoretical values may be used to check the quality of the posterior distributions derived from the model (Deutsch 1997;Olea 2012). These approaches were initially formulated for the estimation/simulation of univariate random functions or for cases where a clear primary variable is to be modeled, with one or more covariates which are of secondary importance. However, in the case of fully multivariate data such as directional (van den Schaeben 2002a, 2002b) or compositional data (van den Boogaart and Tolosana-Delgado 2013; Pawlowsky-Glahn and Egozcue 2020), the entire regionalized vector is seen as an entity which needs to be modeled rather than just its component parts. As a consequence, geostatistical estimation and simulation need to be treated as fully multivariate, and any appraisal of the quality of the geostatistical model needs to take this aspect into account. This concerns all statistical results mentioned before: error statistics, correlation between predictions and observations, and accuracy of estimated intervals.
In this contribution, a generalization of accuracy in the sense of Deutsch (1997) is proposed for the multivariate setting. The proposal is analogous to the method described by Olea (2012) for quantifying the quality of the estimated distribution. Of specific interest is the evaluation of the suitability of a geostatistical estimation or simulation model in the compositional framework, namely, where each variable is non-negative, and its values inform of the relative abundance of a certain component forming the system (Tolosana-Delgado et al. 2019).
After this introduction, four further sections follow. In Sect. 2 the fundamentals of compositional data analysis are recalled briefly along with their implications in geostatistics. Section 3 reviews the existing proposals for univariate validation, with methods, diagrams, and statistics commonly used for this task. In Sect. 4 a fully multivariate approach to cross-validation is proposed for vector-valued random functions specified on the example of compositional data, both for cokriging outcomes (Sect. 4.1) and cosimulation (Sect. 4.2). Two case studies are presented in Sect. 5 illustrating these two sets of techniques. Conclusions are provided in Sect. 6.

Regionalized Compositions and Their Geostatistical Treatment
c, u ∈ A} of compositional data defined on some study region A, where u ∈ A denotes a location in A, and c is an arbitrary but fixed constant. To avoid problems arising out of the fact that compositional data are closed to that constant sum and formed by non-negative components, compositions are usually transformed prior to any geostatistical or statistical treatment. Several logratio transformations are commonly used. These include the centered (clr; Aitchison 1986), additive (alr; Aitchison 1986) and isometric (ilr; Egozcue et al. 2003) logratio transforms. However, the choice of logratio transformation does not impact the final results, because the geostatistical techniques discussed here are affine-equivariant (Filzmoser and Hron 2008;Tolosana-Delgado et al. 2019). Affine equivariance implies that m(ZB) m(Z)B and

S(ZB) B T S(Z)B,
where m denotes the mean, S a covariance matrix, and B a linear transformation. In spite of this invariance, it is mathematically convenient to use the ilr transformation in the geostatistical workflow (Pawlowsky-Glahn and Egozcue 2020). The corresponding regionalized composition of ilr-transformed variables will be denoted by {ζ(u) The image space of the ilr transformation is (D − 1)dimensional Euclidean space. The standard workflow then is known as the principle of working in coordinates: 1. Transform the regionalized composition to logratios using a suitably chosen ilr transformation (Egozcue et al. 2003;Tolosana-Delgado and Mueller 2021) Apply the geostatistical technique to the logratios. 3. Backtransform the geostatistical estimate or realization of the logratio scores, The Aitchison geometry version of the Mahalanobis distance is of fundamental importance for the methods introduced in this paper. With a covariance matrix S as defined above, the (square) Mahalanobis distance between two compositions in the Aitchison geometry is analogous to the (square) Aitchison distance Note that, although d AM is more comfortably defined in terms of the ilr-transformed scores, it is an affine-equivariant quantity, hence intrinsic to the composition and the covariance S, and not dependent on the actual logratio transformation being used. One must merely represent S in the same transformation used for the composition, according to Eq. (1). This is not true of the Aitchison distance (Eq. 3), which holds only in terms of the ilr or clr transformations. With the Aitchison-Mahalanobis distance, one can define the additive logistic normal distribution (AL N ) as the probability model with density proportional to This probability density function is, as required, also affine-equivariant, owing to the determinant being one of the invariants of S.

Cross-Validation, Accuracy and Precision in the Univariate Case
Leave-one-out cross-validation is a well-established tool used to assess a given geostatistical model and to determine the best model from a set of competing models. It is based on kriging, and given a sample data set at each location u, a kriging estimate z * K (u) and corresponding error variance σ 2 * K (u) are derived by removing the location from the set and estimating the value of the variable of interest based on the neighboring data. In jackknifing, a separate validation set is assumed to be available, and the sample data are used to provide estimates or simulated values at the jackknife locations. A more general cross-validation approach, known as n-fold cross-validation, is to partition the data set into n disjoint subsets and then apply jackknifing on each of the subsets based on the remaining sample data. Independently of the validation method actually used, one finally has a paired list of observed and estimated values, as well as their kriging variance.
If the true value at u is z(u), then the error is given by e(u) z(u) − z * K (u), and the squared deviation ratio is defined as .
Averaging these quantities over all sample locations results respectively in a mean error (ME) and mean squared deviation ratio (MSDR). The former and the associated mean square error (MSE) should be close to 0 and the latter close to 1 if the geostatistical estimator is adequately defined (Webster and Oliver 2007). Typical diagnostic plots are a scatterplot of the true values against the estimates, a histogram of the standardized errors and a scatterplot of the standardized errors against the estimates.
Analogously to the linear model, here one expects a tight scatter between estimates and true values, a symmetric histogram of standardized errors with mean close to 0 and variance close to 1 (as a weakened version of standardized normality), and the scatter between standardized errors and estimates showing no correlation (Chilès and Delfiner 2012;Webster and Oliver 2007). Cross-validation is routinely performed in geostatistical practice and used for the appraisal of the variogram and trend model, the local neighborhood and parameters associated with the simulation method applied. Performance measures in addition to ME, MSE and MSDR concern the quality of the local posterior distributions. Their assessment was first discussed in the context of univariate geostatistical simulation by Deutsch (1997) and is based on the coverage of the local distributions. If the local distribution is F u with mean μ(u) and standard deviation σ (u), then the indicator function defined for p ∈ (0, 1] allows measurement of the closeness of the true value to the local mean in the context of a symmetric target distribution. If the simulation algorithm is Gaussian, then the parameters of the local distribution are derived via kriging; otherwise, the local distribution is inferred via the generation of a family of simulated values at the sample location leaving the actual value out. The coverage π ( p) is set equal to the average of i(u, p) over all available locations u.
To determine the accuracy of the model, a further indicator variable a( p) is introduced which is set to 1 if the proportion π ( p) of locations falling into the p-interval exceeds p, and to 0 otherwise. The integral A 1 0 a( p)dp then provides a measure of accuracy. That is, A measures the proportion of exact or over-pessimistic confidence intervals around the estimates. A useful means for appraising the accuracy of the model is a plot of π ( p) against p. An accurate model will result in a plot where the pairs of points ( p, π( p)) fall above the bisector line. Two measures of precision are defined in Deutsch (1997), one restricted to pairs of points ( p, π( p)) falling above the bisector line, called precision and defined as P 1 − 2 1 0 a( p) · (π ( p) − p)dp, and the other, called goodness, given by G 1 − 1 0 (3a( p) − 2)(π ( p) − p)dp. For any value of p for which the point ( p, π( p)) lies below the bisector, the departure of π ( p) from p is penalized in this definition, and high values of G correspond to precise models. Thus, an accurate and precise model has values of A, P and G close to 1. It is important to note that for high accuracy, it suffices that the actual coverage is larger than the nominal one (i.e., p < π( p)), while precision and goodness also reward proximity to the bisector (i.e., | p − π( p)| small). The precision measure P only makes sense in the case of accurate models, while goodness G is more generally useful.
An alternative approach for assessing the quality of the local posterior distributions was introduced by Olea (2012). In his approach, the symmetric p-interval used in Eq. (5) is replaced by the unilateral interval −∞, F −1 u ( p) . For each p ∈ (0, 1) and each u, an indicator variable is defined by putting For each p, the empirical probability p * ( p) is set equal to the average of i O (u, p) over all available locations u and plotted against p. The modeling is optimal if the pairs ( p, p * ( p)) fall on the bisector line, and according to Olea (2012) indicates "perfect global agreement between the modeling of uncertainty and the limited amount of information provided by the sample." In practice, however, there will be deviations from the bisector, and they can be quantified via the maximum absolute deviation and the sum of absolute deviations between the empirical and theoretical probabilities. It should be noted that the functions π (•) and p * (•) are related by

The Compositional Case
Here, the implementation depends on whether or not cokriging is applied to derive the parameters of the local distribution. In this case it is assumed that the composition follows an additive logistic normal distribution as expressed by Eq. (4), at least locally, that is, conditional on the estimated composition. As usual in geostatistics, this assumption cannot be formally tested, owing to the presence of spatial dependence, and is to be considered a modeling choice. Otherwise, if conditional additive logistic normality is deemed inappropriate, one should use either multipoint methods or else some form of multivariate transformation to normality (e.g., Barnett et al. 2014;van den Boogaart et al. 2017;Sepulveda et al., under review) followed by Gaussian cosimulation, in which cases Sect. 4.2 applies.

Cross-Validation, Accuracy and Precision via Cokriging
For compositional data sets, the implementation of the cross-validation procedure via cokriging is straightforward, even if not implemented in most software packages (Tolosana-Delgado and Mueller 2021). At a location to be cross-validated, the entire compositional vector is removed and the surrounding compositions are used to estimate the composition at the sample location via cokriging. As the cokriging is performed in logratio coordinates, the error measures are also calculated in terms of logratios. The mean error is given by and the associated mean square error is MSE , which corresponds to the average Aitchison distance (Eq. 2) between estimates and observations. There are two ways to generalize the mean square deviation ratio to a multivariate quantity (Tolosana-Delgado and Mueller 2021), namely and In the equations above, z(u α ), ζ (u α ), ζ * C K (u α ) and z * C K (u α ) denote the true compositional vector, its logratio image, the logratio estimate and the corresponding backtransform at location u α . The expression ln z(u α )/z * C K (u α ) is an abbreviation of ln z 1 (u α )/z * C K ,1 (u α ) , . . . , ln z D (u α )/z * C K ,D (u α ) . The matrix C K (u α ) denotes the cokriging error variance-covariance matrix at location u α , and σ 2 ii (u α ) are its diagonal elements.
Only the mean error ME is a vectorial quantity. All other quantities (MSE, MSDR 1 and MSDR 2 ) are scalars. The measure MSDR 1 is nothing other than the average over the square Aitchison-Mahalanobis distances (Eq. 2) d 2 AM (z(u α ), z * C K (u α )| C K (u)) between z(u α ) and z * C K (u α ) with respect to the cokriging error variance-covariance matrix. Its target value is equal to D − 1. Moreover, under the hypothesis of additive logistic normality of the D− component compositional random function, the square Aitchison-Mahalanobis distance follows a χ 2 (D − 1) distribution. The version of MSDR in Eq. (7) is the average over the univariate MSDR values for the components of the logratios. In contrast to MSDR 1 , this measure does not have the equivariance property.
The relationship between the Aitchison-Mahalanobis square distances and the χ 2 − distribution gives rise to a diagnostic tool that may be used in place of the histogram of standardized errors of the univariate case. This is a qq-plot of the observed quantiles of the Aitchison-Mahalanobis square distances against the quantiles of the χ 2 − distribution with D − 1 degrees of freedom. As in the univariate case (where the comparison is against the standard normal distribution), one expects the qq-plot to be close to the bisector. Even if multivariate additive logistic normality does not hold locally for the compositional random function, this diagram will still provide a means to rank competing geostatistical parameter setups (mostly variogram models or kriging neighborhoods) in the sense of their approximation to this distributional assumption, exactly in the same way as for the univariate case.
Other common diagnostic plots are scatterplots of the individual components against their estimates, or of the estimation errors versus the estimates. Compositional versions of these plots are formed by the set of scatterplots between pairwise logratios of estimates against pairwise logratios of true values, and a set of scatterplots of logratio estimation errors against logratio estimates (Tolosana-Delgado and Mueller 2021); both can be conveniently presented in a (D × D) matrix of scatterplots.
Compositional analogues of accuracy and goodness described in Sect. 3 are also based on the square Aitchison-Mahalanobis distance of estimates and true values with respect to the estimation error covariance. The definition of the indicator variables required for the calculation of coverage is subject to a certain arbitrariness, because there is no natural ordering of vectorial quantities. In general, any one-dimensional summary of the random composition can be used to generate coverage indicators, as long as the probability distribution of this target quantity is known. A reasonable requirement is for these quantities to be affine-equivariant. The square Mahalanobis distance (Eq. 2) arises naturally as the best option: for each location u and each p ∈ (0, 1], the indicator function i AM (u, p) is defined as analogously to Olea's (2012) proposal. As in the univariate case, the coverage π AM ( p) is defined as the average over all sample locations, and the indicator variable a( p) is equal to 1 for π AM ( p) > p and 0 else. The metrics A, P and G then have the same definitions as previously, and an accurate and precise model has values of A, P and G close to 1. Other one-dimensional summaries can also be of use: It should be noted that the univariate measures discussed in Sect. 3 may be computed for each relevant logratio variable, each one of them being univariate summaries of the random composition. Even the original variables could be considered in this sense appropriate univariate summaries, if one were ready to obtain the confidence intervals in Eq. (5) by means of Hermite quadrature, as explained in Pawlowsky-Glahn and Olea (2004).

Cross-Validation, Accuracy and Precision via Simulation
When cross-validation or jackknifing needs to be based on simulation at the sample locations, the definitions of the errors and also the local distributions are based on the simulation results. As before, the errors are in the first instance calculated in logratio coordinates. For L realizations ζ (u α )| 1, . . . , L, α 1, . . . , N , we define the local mean and the local covariance as Then, analogously to the cokriging case, one has the mean error given as .
In contrast to Eq. (6) where the target value was (D − 1), the value of MSDR 1,sim for a good model should be slightly greater than this quantity, owing to the fact that −1 (u α ) is estimated from a set of L realizations. A reasonable target quantity can be derived from the expected value of MSDR 1,sim under the assumption that ζ(u α ) is normally distributed conditionally on ζ (u α ), which produces a Hotelling's T 2 distribution for the Mahalanobis distance, with parameters (D − 1) and (L − 1). This gives an expected value of thanks to the equivalence between Hotelling's T 2 and Fisher F-distributions and the fact that the expected value of a Fisher F ( p,q) -distributed variate is q/(q − 2). For deriving the coverage indicators to calculate accuracies and derived quantities, we can follow the same ideas as in the section about cokriging: one needs to select a univariate summary of the composition, compute that statistic for the true value and generate its probability distribution with the realizations. Again, in general terms, it makes sense to enforce that summary statistic to be affine-equivariant. The indicator variable defining the position of the true compositional vector within the local distribution can then be based on the square Aitchison-Mahalanobis distance: the distance d 0 (u α ) between the true composition and the mean is compared with the empirical distribution of the distances {d (u α ) : 1, 2, . . . , L} of the simulated results and the mean at u α where (u α ) ζ (u α ) − ζ (u α ) for each 1, 2, . . . , L. The values {d (u α )| 1, . . . , L} are arranged in ascending order d (u α )| 1, . . . , L , and where pL denotes the smallest integer greater than pL. The statistics A, P and G are then defined as in Sect. 4.1. This construction mimics that of Deutsch (1997) in that in the univariate case, the simulated values are used to derive a local distribution, and the location of the true value is determined relative to it. Note that, as in the case of cokriging, other uni-dimensional summaries may also be meaningful for specific cases: the simulation approach outlined here is particularly useful in these cases because the probability distribution of the target summary can always be derived from the set of realizations.

Illustration
In what follows, two applications are provided. The first demonstrates the crossvalidation for cokriging of a regionalized subcomposition of the Tellus data (Young and Donald 2013), while the second concerns cross-validation and accuracy assessment in the case of direct simulation for the modeling of the structure of a tailing storage facility.

The Tellus Data Set
The composition considered in this example (Fig. 1)  The subcomposition was closed through the inclusion of an additional component called Rest and transformed to logratios via the default ilr transform, as described in Tolosana-Delgado et al. (2019). The ilr variables exhibit geometric anisotropy with direction of greatest continuity N135 and a linear model of coregionalization comprised of a nugget, and an exponential structure of range 35 km was fitted, with an anisotropy ratio of 0.4. For the sake of comparison, an isotropic version of this model was also considered, with a range of 26.9 km. Tenfold cross-validation via ordinary cokriging with a moving neighborhood (search radius 60 km, minimum number of samples: 7, maximum number of samples: 20) was applied with both models.
A summary of the performance measures in Table 1 indicates that the anisotropic model is superior to the isotropic one. This is also supported by the graphs of coverage against confidence level in Fig. 2. Additionally, the actual coverage of the compositional model is greater than theoretically expected for the majority of the confidence levels, as evidenced by the value of A compared to those of the individual components. The accuracy plots in Fig. 2 provide further insight on the behavior of coverage versus confidence for the components and the entire composition. The coverage is generally closer to the chosen confidence level for the individual components compared to the model in its entirety (although the overall accuracy of the compositional model is superior to those of the constituent parts) for values of confidence ( p) up to 0.7; for greater values of confidence, this behavior is no longer observed. This example thus provides support for our claim that evaluation of the performance of a geostatistical model for compositional data should not be based on the performance of the model on the individual components only. Figure 2 also shows that the model ignoring anisotropy appears to have higher accuracy, but at the price of constructing much wider p-intervals. In the accuracy plot, this is shown through the much greater deviation from the bisector. Accounting for anisotropy notably increases the model precision, as evidenced by the higher value of P 0.79 and G 0.88 in the anisotropic case compared to P 0.47 and G 0.74 for the isotropic model.   . 2 Accuracy plot derived from tenfold ordinary cokriging cross-validation of the Tellus subcomposition sample, with an isotropic model (omni) and with an anistropic one (globally, and for each one of the ilr variables based on the anisotropic LMC) Table 2 gives a summary of the compositional error measures for both the isotropic and anisotropic models. The vector of ilr mean errors (ME) is in both cases very close to zero, with mean square errors (MSE) of roughly 0.4 in both cases. The difference between the two models is apparent in the MSDR measures, where it is clear that the isotropic model overestimates the spread. MSDR 2 clearly shows that on average over all logratios, the anisotropic model meets the target of MSDR equal to 1.
This subcomposition was also studied using direct simulation (Mariethoz et al. 2010), with a validation subset of 100 samples, and the rest used as conditioning data and training image generation, reported in Tolosana-Delgado and Mueller (2021; Ch. 10-11). To complement the discussion in these chapters, the accuracy metrics in the original units (raw data) are shown in Table 2. The results show overly optimistic metrics in general terms, indicating a potential lack of variability in the training image,  (Table 3).

The Tailings Storage Facility Case Study
The second example illustrates the usage of the simulation-based approach of Sect. 4.2, by means of a case study with direct sampling (Mariethoz et al. 2010) with a very complex setup, far from additive logistic normality. Here the aim is to model distributions of particle types in a multiple stream tailings storage facility (Selia et al.  in prep). Stratigraphic forward modeling is used for training image generation and direct sampling for data fusion with the measured data. A multipoint method was considered as most appropriate because of the strong effects of non-linearity and nonstationarity of the patterns typical of these anthropogenic sedimentary systems. Briefly, the study considers four particle classes, according to size and dominant mineralogy, forming a four-part mass composition: sand-sized silicates (V1), clay-sized silicates (V2), sand-sized sulfides (V3) and clay-sized sulfides (V4). The forward model used is Delft3D-FLOW (Lesser et al. 2004), an open-source, process-based stratigraphic forward modeling software accounting for diffusion, advection in both bed load and  suspended load, erosion and compaction in both aerial and subaquatic environments. The parameters of the forward simulation are described in Table 4 (left). Boundary conditions are designed to mimic the behavior of tailings dams, allowing water to seep while retaining sediments. Results were cropped to the basin without the upstream channels and upscaled to 18 × 21 × 21 voxels to form the ground truth, which for this study, is also taken to be the training image (Fig. 3). Synthetic boreholes were randomly taken at 36 locations (Fig. 4). Of the resulting 677 samples, 100 were randomly picked for leave-one-out cross-validation.
To predict each of these 100 samples, direct sampling was applied excluding that sample, based on the parameters specified in Table 4 (right), using the Aitchison distance (Eq. 3) as the measure of proximity between the composition at each training image pixel and the data set. The resulting simulations were used to calculate the accuracy and goodness as defined in Sect. 4.2. As can be seen from the numerical (Table 5) and graphical results (Fig. 5), the simulations show inaccurate but moderately good results: the accuracy curve for both individual ilr coefficients and for the global Aitchison-Mahalanobis summary are systematically below but close to the reference line, particularly for theoretical coverage levels below 0.20. This indicates that confidence intervals tend to be too small to deliver the coverage promised by their nominal confidence. Correspondingly, the numerical values of accuracy A are close to zero, and precision is meaningless (hence not reported in the tables). In this situation, the goodness G becomes useful: G is above 0.81 for all three ilr variables, and the global goodness is the highest (0.89).
The availability of a simulation allows an evaluation of the accuracy in terms of the original four components. Results are reported in Table 6 and Fig. 6, and do not qualitatively diverge from the logratio-based results: confidence intervals are too large, so that accuracy values A are not useful. However, goodness values are high (above 0.83), suggesting that the model is acceptable in both representations, raw and logratio.

Conclusions
We have provided an extension of quantitative measures of goodness of fit to the multivariate context, particularly for compositional data, which allow ranking of models analogously to already existing univariate approaches. In addition to generalizations of mean errors to vectorial quantities, the average Mahalanobis distance between estimates and observations (MSDR 1 ) gives an additional decision tool for choosing between competing models. This measure explicitly accounts for the multivariate structure, is affine-equivariant and shifts the focus from the individual components to the entire composition. Alternatively, an average mean square deviation ratio over all logratios (MSDR 2 ) can also be used, although this quantity is not affine-equivariant, that is, it depends on which specific logratios are being used to compute the statistic. In principle, both measures focus on different aspects and could result in different rankings of models, although both would result in the choice of the same model in the case studies provided here. The joint accuracy and precision measures introduced here can provide further insights into the system inasmuch as they evaluate the global structure, and may rank competing models differently as when evaluated in terms of just the marginals.
The Mahalanobis distance also provides a reasonable one-dimensional summary of the multivariate distribution, to derive a cumulative distribution and with them measures of accuracy, precision and goodness after Deutsch (1997). Accuracy cannot be evaluated isolated from goodness or precision measures, in either the univariate or multivariate case. This is particularly evident in the metrics for the isotropic versus anisotropic linear models of coregionalization of the Tellus subcomposition, where the accuracy is marginally higher for the isotropic model, but goodness and precision provide a higher contrast to choose the better model. In particular, the goodness metric is important here, as it accounts for both accurate and inaccurate cases. In the case of the tailings storage facility case study, the model overall ended up showing results being too far from the center and hence resulting in an accuracy value of 0. Nevertheless, examination of the actual versus theoretical coverage shows reasonable goodness, with the plot for the global metric quite close to the bisector. The tools could thus still be used for ranking competing models, such as training images or method parameter setups, had several of them been available.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.