## Abstract

Traditional uses of data transformations in geoscience have typically been motivated by three objectives: (1) creating normally distributed data; (2) creating data that are additive; and (3) making errors constant across the data range. Part 1 of this paper discussed the transformation applicability for the first two reasons in geochemical applications, and introduced a fourth motivation whereby transformation maximizes the variance (or geochemical contrast) in the transformed data; it did not discuss transformation to stabilize measurement error.

Although transformation to stabilizing errors in geochemical data is not common, this is a useful attribute in geochemical data analysis. The transformation is dependent on the model describing the magnitude of measurement (sampling and analysis) error as a function of concentration. Models describing measurement error as a function of concentration can be used to derive a transformation that will stabilize the measurement error in the transformed variable. Poisson, binomial and hypergeometric models are typically used to describe sampling errors, whereas straight line (constant, proportional and affine) models are used to describe analytical errors. The associated variance stabilizing transformations, derived from these models, have constant propagated errors. As a result, these transformations create a ‘level playing field’ for subsequent data analysis, enabling the discovery of additional information in the data.

Homoscedastic measurement error allows the geochemist to justify use of a specific transformation based not on the subsequent data analysis results (circular reasoning; ‘the end justifies the means’), but on optimal properties created by the transformation. In this way, objective results can be achieved scientifically, providing another motivation to collect geochemical quality control data.

- geochemistry
- data analysis
- transformation
- constant variance
- Poisson
- binomial
- hypergeometric
- measurement error
- distribution

## Introduction

In Part 1 of this paper, Stanley (2006) reviewed the three traditional motivations for transformation of geochemical data – (1) normality, (2) additivity and (3) homoscedasticity – and introduced a fourth motivation whereby a power transformation is used to maximize data variance (geochemical contrast). Stanley (2006) indicated that transformation for the purpose of creating normally distributed data is probably misguided, because most geochemical datasets are multi-modal, and thus cannot be ‘normalized’ with a simple monotonic transform that does not fundamentally alter the nature of the data. Similarly, Stanley pointed out that additivity is a characteristic that commonly exists in trace element geochemical concentrations, and so transformation to create additivity is seldom necessary unless major oxide concentrations are being considered. However, in his paper, Stanley (2006) omitted detailed discussion of transformation for the purpose of creating homoscedasticity, a feature whereby measurement error is constant across the range of the transformed variable. As a result, a discussion of the nature of measurement error and how it changes across the range of element concentrations in a dataset, as well as how one can identify a transformation that will propagate this heteroscedastic error in the original concentration data into homoscedastic error in the transformed result, is presented below.

## Error variance stabilization

‘Error variance stabilization’, or the creation of homoscedastic error, historically has not been a traditional motivation for data transformation in geochemical applications, but is a desirable trait for geochemical data that are to be evaluated numerically or statistically. Many statistical procedures assume equal amounts of error in each observation (e.g. ordinary least squares regression). As a result, error variance stabilization transformation can be used to treat data that are ‘heteroscedastic’, having magnitudes of variations (measurement errors) that change across the range of the data.

Most geochemical concentration data are ‘heteroscedastic’, as sampling and analytical (or measurement) errors commonly increase with concentration (Thompson & Howarth 1973, 1976*a*, *b*, 1978; Thompson 1973, 1982; Fletcher 1981; Stanley & Sinclair 1986; Stanley 2003*b*). This imposes a complication to data evaluation, because the variations attributable to measurement error change their magnitude across the range of the data, and the criteria used to determine whether observed differences are real thus change across the range of concentrations as well.

Transformation to stabilize the error in a geochemical variable creates a new variable that is ‘homoscedastic’ (has constant error variance across the range of the variable), and this attribute can significantly simplify subsequent statistical tests, numerical procedures and the graphical presentation of the data. The transformation function employed must be specific to the model that describes how error changes with concentration. Propagation of the heteroscedastic error described by this model through the associated transformation into the new variable results in new error that is homoscedastic. Subsequent evaluation of the new homoscedastic data is likely to be far simpler and may lead to more compelling conclusions because the data can be evaluated on an ‘even playing field’.

Examples of error stabilizing transformations commonly used in other scientific applications (Hoyle 1973) include (1) the logarithmic transform (for proportional error): (1) (2) the square root transform (for Poisson error): (2) and (3) the angular transform (for binomial error): (3)

To illustrate how an error stabilizing transformation works, consider case 2. If Poisson measurement error exists, then *σ _{c}*= (from the variance of a Poisson distribution:

*σ*

^{2}

_{x}=

*x*; Speigel 1975), where

*c*is an element concentration. Propagation of this error into the square root of

*c*(

*y*=) requires derivation of the associated error propagation formula, which itself requires use of the generative error propagation equation (Stanley 1990): (4) In this case, there is only one variable (

*c*), which we will assign to

*x*(

_{1}*p*=1), so Equation 4 simplifies significantly. The only required partial derivative is: (5) and thus Equation 4 becomes: (6) As a result, regardless of the magnitude of

*c*, the propagated measurement error in

*y*(the transformed variable) equals ¼ and is constant (homoscedastic). Stanley (2006) demonstrated that proportional measurement errors become homoscedastic when data are transformed using a logarithmic function, and analogous procedures can be used to demonstrate that binomial errors become homoscedastic when data are transformed using an angular function.

To date, no published geochemistry-related paper is known to the author that transforms geochemical data for this purpose, although ‘homoscedasticity’ would significantly facilitate a significant number of efforts to evaluate geochemical data.

### Transformations to stabilize measurement error

In geochemistry, compositional variation can be introduced by a number of sources. The variation of interest is commonly introduced by geological processes (*g*); however, additional variation is also commonly introduced during procedures used to determine the compositions of geological materials: sampling (*s*) and geochemical analysis (*a*). Thus, the total concentration variance (*t*) in a suite of geological samples is the sum of these three variances: (7) where *c* refers to the element concentrations. The latter two variance terms are commonly not of ultimate interest to a geochemist, and are thus here referred to collectively as measurement errors (*mc*): (8)

Geochemical measurement errors interfere with our ability to understand the underlying geochemical variations (*σ _{gc}*

^{2}) present in a suite of geological samples because they impose additional variation that obscures the compositional variations caused by geological processes. Unfortunately, these error variances typically affect the true element concentrations in a suite of rocks differently across the range of concentrations because the magnitudes of these errors are commonly a function of concentration. Thus, larger concentrations often have larger measurement errors, and

*vice versa*, and the magnitude of measurement error changes across the range of the data.

This heteroscedastic variance attribute has important implications in subsequent data analysis. Consider two samples with Pb concentrations of 5 and 10 ppm. Because these relatively low concentrations would typically have low measurement errors (albeit high relative errors), these two concentrations may be significantly different ( Fig. 1), a feature that could be demonstrated using a statistical test (e.g. a *t*-test). However, this may not be true for two samples with higher Pb concentrations that exhibit the same concentration difference, say 35 and 40 ppm. This is because larger measurement errors would typically be associated with those analyses, as the two concentrations may not be statistically different given measurement error. As a result, samples with low concentrations have to be evaluated differently from samples with high concentrations, and this adds a complication to data evaluation.

### Sampling errors

Because geochemical data have non-constant measurement errors, and thus typically do not lie on an ‘even playing field’, it is worth considering how the magnitude of geochemical error changes with concentration. The error–measurement relationship can be modelled in two different ways: theoretically or empirically. Sampling error has historically been described theoretically because a simple ‘equant grain model’ can be invoked to explain sampling error in most samples (Clifton *et al*. 1969; Ingamells & Switzer 1973; Ingamells 1974*a*, *b*, 1981; Cheng 1995; Stanley 1998, 2003*a*). In the ‘equant grain model’, mineral grains are thought to be of equal size and shape, and to consist of either an ore mineral or gangue mineral of constant composition. With this simple model, the sampling variance structure can be modelled ideally, because simple statistical distributions can be used to describe the sampling errors in real samples. Three statistical distributions have traditionally been employed to estimate the magnitude of sampling errors using the ‘equant grain model’: binomial (Wickman 1962; Kleeman 1967; Visman 1969; Cheng 1995;Stanley 2003*a*), hypergeometric (Stanley 2003*a*) and Poisson (Clifton *et al*. 1969; Ingamells & Switzer 1973; Ingamells 1974*a*, *b*,Stanley 1984) distributions.

Real samples clearly do not have mineral grains of ideal size, shape and composition. As a result, they exhibit a more complex sampling variance structure. The ‘equant grain model’ can be used to describe this real sampling variance structure using ‘effective’ (ideal) sampling parameters. The magnitudes of these ‘effective’ sampling parameters can be determined through regression against real replicate sample results. The regressions provide estimates of the ‘effective’ sampling parameters that define a sampling variance structure for the ‘equant grain model’ that replicates the sampling variance structure in the real samples. As a result, in spite of the overall mineralogical complexity observed in the real samples, ‘effective’ sampling parameters can be used to quantitatively describe the magnitude of the real sampling error as a function of concentration. The character of these ‘effective’ sampling parameters depends on the statistical distribution employed, which in turn depends on the nature and geochemical behaviour of the element of interest in the samples.

#### Poisson sampling errors

For example, if we consider samples containing a small number of nuggets (typical of gold ores), the Poisson distribution (*p*) is the appropriate distribution to be employed (Garwood 1936; Clifton *et al*. 1969; Ingamells & Switzer 1973; Ingamells 1974*a*, *b*, 1984; Stanley 1998, 2003*a*). This distribution applies to cases where rare events occur over time or space (Speigel 1975). As a result, the collection of a sample with a small number of gold nuggets is a case where the Poisson distribution effectively describes the sampling error.

Unlike the normal distribution, the Poisson distribution has a variance that is functionally related to the observed number of rare events (nuggets): (9) where *x* is the number of rare grains collected in a sample, and *σ _{px}*

^{2}is the Poisson sampling variance for this rare number of grains.

Unfortunately, geochemical analyses do not typically report the number of nuggets in a sample, but rather report the concentration of an element in a sample, as a result of inclusion of those nuggets (e.g. the gold concentration). Thus, Equation 9 cannot be used to describe the sampling error of element concentrations in nugget-bearing samples, because it does not relate sampling error to concentration. However, if one manipulates Equation 9 by inverting it, and then multiplying both sides by the standard deviation (*σ _{px}*), we obtain: (10) The relative sampling error on the number of nuggets in a sample (

*σ*/

_{px}*x*) must equal the relative sampling error on the element concentration (

*σ*/

_{pc}*c*; both parameters are unitless; Clifton

*et al*. 1969; Stanley 1998, 2003

*a*): (11) Hence, the relative error of replicate concentrations of the nugget-borne element can be used to determine the ‘effective number of nuggets’ in a sample (Clifton

*et al*. 1969; Stanley 1998, 2003

*a*).

Once the ‘effective number of nuggets’ (*x*) in a set of replicate samples has been determined, the average concentration of the nugget-borne element of interest in those samples (c̄), the mass of the replicate samples (*M _{s}*), the average nugget shape in those samples (here assumed to be spherical), and the nugget density (

*ρ*), along with the ‘effective number of nuggets’ (

*x*, from Equation 10) can be used to determine the ‘effective size of the nuggets’ (here expressed as a diameter,

*d*; Clifton

_{x}*et al*. 1969; Stanley 1998): (12)

Furthermore, the average concentration (c̄) and the ‘effective number of nuggets’ (*x*) can be used to determine the ‘effective concentration contribution per nugget’ (*k*): (13)

The ‘effective concentration contribution per nugget’ (*k*), like the ‘effective size of the nuggets’ (*d _{x}*), is also a fundamental ‘effective’ sampling parameter that establishes the relationship between the ‘effective number of nuggets’ in a sample and the sample concentration. It likewise establishes the relationship between the sampling standard deviation on the ‘effectivenumber of nuggets’ and the sampling standard deviation on the concentration: (14)

Substituting the relationships described in Equation 14 into Equation 9 and rearranging, we obtain the relationship between the concentration sampling error, produced by the random collection of nuggets in a sample, and the resulting concentration: (15) Thus, the ‘effective concentration contribution per nugget’ value (*k*) is a fundamental sampling parameter that can be used to quantitatively describe a Poisson sampling error structure (the measurement–error relationship) that matches the observed error structure observed in the replicate samples ( Fig. 2; Stanley 1998, 2003*a*). Clearly, in samples containing rare nuggets, the associated Poisson sampling error defines a heteroscedastic relationship where error varies across the range of data (an ‘uneven playing field’; according to Equation 15) on which the geochemist is forced to undertake data interpretation for elements contained in rare nuggets.

#### Binomial sampling errors

If an element occurs in a mineral grain that is fairly common in a suite of samples, the corresponding sampling error structure cannot be described by a Poisson distribution. Rather, an alternative distribution must be considered. Because an element concentration in a sample is technically equal to the probability of randomly selecting an atom of the element of interest from the sample (if this were physically possible; Stanley 2003*a*), there are effectively two types of atoms in a sample: atoms of the element of interest (Cu, say), and atoms of elements that are not of interest (‘non-Cu’). Sampling (collecting) a large number of these atoms at random to produce a ‘sample’, provided that there is a sufficiently large population of atoms from which to collect the sample so that random selection of each atom without replacement does not alter the selection probabilities, would result in sampling errors that are binomially distributed (*b*; Wickman 1962; Kleeman 1967; Visman 1969; Speigel 1975; Cheng 1995; Stanley 2003*a*; note that ‘binomial’ means ‘two names’, in this case ‘Cu’ and ‘non-Cu’).

In most geological materials, elements occur in mineral grains. Just like samples of atoms in the above example, samples of mineral grains will also be subject to binomially distributed sampling errors (e.g. ‘chalcopyrite’ and ‘non-chalcopyrite’), at least provided that (i) the equant grain model applies, and (ii) samples are sufficiently small (in mass) relative to the populations from which they are drawn (Speigel 1975). The magnitudes of these sampling errors for each mineral grain type can be predicted using the formula for the variance of a binomial distribution (Stanley 2003*a*): (16) where *x* is the number of mineral grains containing the element of interest (‘chalcopyrite’), *n–x* is the number of mineral grains that do not contain the element of interest (gangue; ‘non-chalcopyrite’), *c* is the concentration (expressed as a proportion) of the mineral grain containing the element of interest (chalcopyrite; *c=x/n*), and *n* is the total number of mineral grains in the sample ( Fig. 3).

Unfortunately, elements may reside in more than one mineral grain type within a sample (e.g. chalcopyrite, bornite and chalcocite) that may differ in size and composition. As a result, the sampling errors on the element concentrations may not be truly binomially distributed; rather they may exhibit a distribution resulting from samples that are physical mixtures of several mineral grain types of different sizes and compositions containing the element of interest, where each of the mineral grain types has sampling errors that may be binomially distributed. Although this underlying element concentration sampling error distribution is unknown, one can still use the ‘equant grain model’ to describe an ideal sample that contains the same variance structure as the real samples. Employing the ‘equant grain model’ for common mineral grains is analogous to the approach used to describe Poisson sampling error in samples containing rare mineral grains, above (Clifton *et al*. 1969; Ingamells 1974*a*, *b*; Stanley 1998, 2003*a*).

Because elements typically occur in several types of mineral grains within a sample, the value of *n* obtained by applying the ‘equant grain model’ to the variance structure of samples with common mineral grains of interest, will not be the actual number of mineral grains in the sample, but rather will describe the ‘effective total number of grains’ in the sample. In other words, the relationship between the observed sampling variance and concentration (defined by *n*) in the real samples is equivalent to that of an equant grain model with a certain number of ‘effective’ grains (*x*) containing the element of interest, and with *n* ‘effective’ total number of grains). Consequently, *n*, like *k* for the Poisson sampling error, is a sampling parameter, in this case for binomial sampling error, and can be used to describe the relative variability in sampling due to binomial sampling errors of elements contained within common mineral grains in a sample.

The non-linear relationship defined by Equation 16 illustrates that element concentrations in common minerals contained in samples will exhibit heteroscedastic errors, and thus will force the geochemist to evaluate the concentration data on a tilted playing field. This binomial variance function also guarantees that sampling errors will be largest at 50% concentration and lowest (=0) at 0 and 100% concentration (Fig. 3). Furthermore, the larger the value of *n*(the ‘effective total number of mineral grains’), the smaller the binomial sampling error will be, and the closer to normal these binomial sampling errors will be distributed (Speigel 1975). This is because, as *n* → ∞, the binomial distribution converges on the normal distribution. Nevertheless, even with large *n*, binomial sampling errors still retain a dependency on the concentration level of the sample according to Equation 16, unlike normal sampling errors, which have no such constraint.

#### Hypergeometric sampling error

In the above discussion, it was assumed that the size of the sample of mineral grains was small relative to the size of the population from which it was collected. This was necessary to guarantee that the probability of randomly selecting a mineral grain of interest in a sample does not change in the parent material even if the grain is not returned to the sample after selection (geological samples are not returned to the sample site, i.e. they are not replaced; Stanley 2003*a*). Effectively, collection of a sample composed of many grains is a not a binomial process, because after selection, each grain is not replaced. However, when the size of the population from which a sample is drawn is large, relative to the size of the sample (by at least two orders of magnitude), then even if the selected mineral grain is not replaced, the probability of selecting a second grain containing the element of interest (and adding it to a sample) does not change significantly (Stanley 2003*a*), and can be described by a binomial distribution.

Unfortunately, if the size of the population is small relative to the size of the sample collected from it, the probability will change, and so an alternative distribution must be considered. The distribution that applies to cases involving ‘sampling without replacement’ is the hypergeometric distribution (*h*; Speigel 1975; Stanley 2003*a*). This distribution is fundamentally related to the binomial distribution, and thus has a variance that is similar to that of a binomial distribution: (17) where *n* and *c* are as before, and *N* is the size of the population ( Fig. 4). When the size of the population is very large relative to the size of the sample, then (*N* − *n*) and (*N*− 1) are approximately equal, and the hypergeometric variance simplifies to the binomial variance (Equation 16; Stanley 2003*a*).

Unfortunately, in the formulation of Equation 17, there are two sampling parameters (*N* and *n*), making it difficult to describe the sampling error structure uniquely (using a single sampling parameter). As a result, Equation 17 can be modified by dividing the numerator and denominator of the extra ‘hypergeometric’ term [(*N*− *n*)/(*N*− 1)] by *n*. This yields: (18) where *d*=*N*/*n*.

In this model, *d* is a parameter, because it is merely the ratio of the size of the population from which a sample is collected, divided by the size of the sample. Thus, *d* is effectively a measure of the inverse of the sampling density of a geochemical survey, and can be determined from the survey grid geometry and sample size. For example, if soil samples on a 50-m by 50-m grid have a volume equivalent to a cube 10 cm on a side (corresponding to a sample of approximately 2 kg in mass), and if we assume that the appropriate soil horizon from which the samples could be collected is continuous across the 50-m by 50-m sampling area from which each sample is collected, and is 10 cm thick, the size of the population from which the sample could be collected is 50 m by 50 m by 10 cm (or 250 m^{3}). Because the sample is 0.001 m^{3}, *N*/*n*=*d*=250 000. In this case, *d* is large enough (>100) so that the hypergeometric and binomial variances will be nearly identical.

However, if geochemical samples from drill core consist of 50-cm pieces of one-half of NQ core (corresponding to samples of *c.* 4 kg mass), and are collected at intervals of once every 10 m, *d* is much smaller (=40), and the hypergeometric variance will be very much smaller than the binomial variance.

Utilizing the ‘survey density parameter’ *d* in the hypergeometric variance formula causes Equation 18, like Equation 16, to have only one sampling parameter (*n*), and thus, like the Poisson and binomial models, an equant grain hypergeometric model can be employed analogously to describe the sampling variance structure.

Note that the Poisson, binomial and hypergeometric variance formulae of Equations 15, 16 and 18 are fundamentally different from that of a normal distribution. This is because a normal distribution does not require that there be a functional relationship between the magnitude of the mean and the magnitude of the standard deviation (for a normal distribution with a given mean, the standard deviation is unconstrained). Instead, all three sampling models impose a functional relationship between concentration and variance (Equations 15, 16 and 18) that ensures that the sampling error is not constant across the range of the data.

Furthermore, by examining all three Poisson, binomial and hypergeometric models for sampling error (Equations 15, 16 and 18), one can demonstrate that each model has an identical variance structure of the form *σ*^{2}=*r*×*f*(*c*), where *r* is a function of the appropriate sampling parameter. Regression can thus be used to estimate the appropriate sampling parameter if replicate measurements are available to determine the concentrations [*c* and thus *f*(*c*)], and the errors () in a suite of samples. Estimation of *r* is achieved by calculating the ordinary least squares estimate of *r*=Σ*xy*/Σ*x*^{2}, and this allows the appropriate sampling parameter (*k* or *n*) to be estimated directly.

### Analytical errors

In contrast to sampling errors, errors introduced during geochemical analysis probably have a more complex relationship with concentration. Analytical procedures introduce errors from a variety of sources (e.g. counting errors, calibration errors, regression errors, etc.), and it is difficult to predict or quantify the actual amount of error from each source to allow a rigorous treatment of the errors. Because most geologists do not have the resources to precisely establish the complex analytical error relationships with concentration, the error–measurement relationship is commonly estimated empirically. Generally, a straight line relationship between concentration and the analytical error standard deviation is assumed, at least in the absence of specific information to the contrary.

A common technique for describing the magnitude of analytical error as a function of concentration involves the use of Thompson–Howarth plots (Thompson & Howarth 1973; 1976*a*, *b*1978, Thompson 1973, 1982; Fletcher 1981; Stanley & Sinclair 1986; Stanley 2003*b*): (19) where *m* is the relative error term (slope) and *b* is the absolute (or constant) error term (intercept). In many cases, only proportional (*m*) error or only fixed (*b*) error is assumed. Equation 19 accommodates the more general case where both relative and absolute errors contribute to the overall analytical error. As with the sampling error models, this ‘affine’ analytical error model also varies with concentration, and thus establishes an ‘uneven playing field’ (i.e. the error is heteroscedastic; note that the term ‘affine’ is used here to describe the general case because models involving proportional, fixed and proportional and fixed error terms are all ‘linear’).

### Measurement–error models

Total measurement–error models for geochemical concentrations are far more complex than the above analytical or sampling error models because they are the sum of the variances associated with these simple models. As a result, the sampling and analytical models can be combined to produce a variety of total measurement–error–concentration relationships. Total measurement–error models describing the various combinations of binomial, hypergeometric and Poisson sampling errors, and proportional, fixed and affine analytical errors are presented in Table 1. All of these measurement–error models are functions of concentration (except the fixed error model, which is already homoscedastic), and thus have errors that change across the range of the concentrations.

To avoid the difficulties that result from having non-constant sampling and analytical errors, data subject to heteroscedastic error can be transformed into a new variable in which the effects of measurement error, propagated through the transformation, are constant across the range of the data (i.e. are homoscedastic). This may be done using error variance stabilizing transformations (Hoyle 1973).

## Error variance stabilizing transformations

Error variance stabilizing transformation functions *g(x)* satisfy the calculus chain rule property that: (20) where *z* is assigned a value of 1, and is the functional relationship between a variable (an element concentration) and its error variance (Hoyle 1973). Equation 20 is identical in form to the formula for error propagation through a function of one variable (a simplification of Equation 4; Stanley 1990): (21) and is set to 1. Consequently, taking the square root of both sides of Equation 20, dividing by the variance of *c*, and integrating with respect to *c* yields: (22)

This equation may be used to determine the transformation that will stabilize the measurement error described by the functional form of (e.g. the total measurement–error models in Table 1). The function *g*(*c*) is thus a transformation that causes the error in *c* to propagate into *g*(*c*) such that the new error equals unity. For example, substituting a total (combined) error variance model for binomial sampling and proportional analytical error (Table 1, Equation T10, column 1) into Equation 22 yields: (23) Solving this integral gives:if *m*^{2}*n*>1 if *m*^{2}*n*<1 and if *m*^{2}*n*=1 (24) for all non-negative *m* and *n* (Beyer 1981). Propagation of the error in *c* (*σ _{mc}*) into the function

*g*(

*c*) will result in an error in

*g*(

*c*) that is constant and equal to one.

The resulting variance stabilizing transformation (Equation 24), for combined binomial sampling and proportional analytical error, is dependent on the magnitude of *m*^{2}*n*, which is essentially a measure of the relative contribution of sampling and analytical error. When *m*^{2}*n*>1, the relatively large values of *m* and *n* describe large analytical and small sampling errors, so analytical error is dominant; conversely, when *m*^{2}*n*<1, the relatively small values of *m* and *n* describe small analytical and large sampling errors, so sampling error is dominant. When sampling error is dominant (*m*^{2}*n*<1), the transformation involves an arcsine function, and thus is similar in form to the variance stabilizing transformation for binomial sampling error alone (Table 1, Equation 3 and T4). Alternatively, when analytical error is dominant (*m*^{2}*n*>1), the transformation involves a logarithmic function, and thus is similar in form to the variance stabilizing transformation for proportional analytical error (Table 1, Equation 1 and T2).

Furthermore, as *m* goes to zero (causing analytical error to diminish to zero), the transformation converges on the angular transform, the error variance stabilizing transform for binomial sampling error (Table 1, Equation 3 and T4,), whereas as *n* goes to infinity (causing the sampling error to diminish to zero), the transformation converges on the error variance stabilizing transformation for proportional analytical error (Table 1, Equation 1 and T2; R. Karsten, 2004, pers. comm.).

As a result, this combined variance stabilizing function (Equation 24) has the capability of levelling the combined effects of both binomial sampling and proportional analytical error in a gradational manner for the continuum of cases involving only analytical error, to cases involving both sampling and analytical error, to cases involving only sampling error. The actual form of the transform is dependent on the relative magnitudes of the two error parameters (*m* and *n*), and thus can be tailored (tuned) to suit different geochemical datasets with different sampling and analytical error magnitudes.

### Estimating sampling and analytical errors

The *n*, *k*, *m* and *b* parameters for the error variance stabilizing transformations under consideration may be estimated for a suite of samples using sample and site replicates. Replicate analyses of splits of crushed and homogenized samples exhibiting a range of concentrations can be used to establish the functional relationship between concentration and analytical error (a measurement–error model; Thompson & Howarth 1973, 1976*a*, *b*, 1978; Thompson 1973, 1882; Fletcher 1981; Stanley & Sinclair 1986; Stanley 2003*b*,⇓). The *m* and *b* parameters are the slope and intercept of the regression line obtained when the mean concentrations of analytical replicates are plotted against the standard deviations of the analytical replicates using the modified Thompson–Howarth technique (Stanley 2003*b*). These parameters allow estimation of the magnitude of analytical error at any concentration.

Similarly, geochemical analyses of replicate samples collected at several sites can be used to estimate the total measurement error. By subtracting the associated analytical error variance from the total error variance for each site replicate, the sampling error variance can be determined at any concentration. The sampling error magnitudes for each site replicate can then be used in Poisson, binomial and/or hypergeometric sampling error models to obtain functions describing sampling error as a function of concentration.

For Poisson error (Equation 15) a least squares regression can be used to obtain a maximum likelihood estimate of the unknown parameter *k* for the site replicate samples. The resulting value of *k* serves as a Poisson sampling parameter for the suite of samples under consideration, and describes the ‘concentration contribution per nugget’ in the equant grain model (Clifton *et al*. 1969; Ingamells 1974*a*, *b*; Stanley 1998, 2003*a*). This ideal Poisson distribution-based equant grain model with nuggets that contribute a concentration *k* per nugget will exhibit the same sampling error at a given concentration as the replicate data under consideration.

For binomial error (Equation 16) a least squares regression can also be used to obtain a maximum likelihood estimate of the unknown parameter *n* for the site replicate samples. The parameter *n* serves as a binomial sampling parameter for the suite of samples under consideration, describing the number of equant grains in samples from an equant grain model (Stanley 2003*a*). As above, this ideal binomial distribution-based equant grain model with *n* grains per sample will exhibit the same sampling error at a given concentration as the replicate data under consideration.

Finally, for hypergeometric error (Equation 18) a least squares regression can also be used to obtain a maximum likelihood estimate of *r*=(*d*−1)/(*nd*−1) for the site replicate samples. Using an appropriate value of *d* derived from the survey grid spacing and sample size (Stanley 2003*a*) allows subsequent calculation of *n*, which serves as a hypergeometric sampling parameter for the suite of samples under consideration (Stanley 2003*a*). As above, this ideal hypergeometric distribution-based equant grain model with *n* grains per sample will exhibit the same sampling error at a given concentration as the replicate data under consideration.

The regression estimates of *n* and *k* from the above regressions of site replicate samples can be used to estimate the sampling error for all samples. Along with the estimates of *b* and *m* obtained from Thompson–Howarth analysis that describe a straight line model for the analytical error, these sampling error models can be used in an appropriate error variance stabilizing transformation (Table 1), to convert the concentrations under consideration into a new variable that has constant measurement error. By using the estimates of *m*, *b*, *n* and *k* derived from a routine quality control evaluation of the geochemical concentration data, the error variance stabilizing transformation can be ‘tuned’ to accommodate the right proportion of sampling and analytical error in the samples.

Below, several example data sets are used to illustrate this procedure.

### Example transformations

#### Poisson sampling error model

An example of the derivation and use of a variance stabilizing transform is presented below using quadruplicate reverse circulation drill core samples analysed for Au from an anonymous Colorado gold prospect. The drill hole database consists of 40 drill holes representing over 25 000 feet (7620 m) in length. More than 5000 samples, 5 feet (1.524 m) in length, were collected, and the first 2500 drill hole samples were analysed in quadruplicate because of the suspicion of a significant nugget effect (the remaining samples were analysed in duplicate, after recognition that a significant nugget effect was not present). Analysis was by classical fire assay using a 1 assay tonne (29.166 g) sample mass. Sample grades were reported in ounces per ton (opt) with a detection limit of 0.001 opt.

The database investigated includes approximately every tenth sample (476 samples) analysed in quadruplicate, plus all quadruplicate samples with an average Au grade >0.02 opt (676 samples) for a total of 1152 quadruplicate samples ( Fig. 5). Unfortunately, no independent assessment of actual analytical error is available that would allow consideration of both sampling and analytical components in the measurement–error model. As a result, the resulting function describing measurement error is assumed to involve only Poisson sampling error, as, in the absence of analytical replicate results, the magnitude of analytical error is assumed to be very small relative to sampling error.

Regressing the Poisson error model (Equation 15) against the replicate data, a maximum likelihood value for *k* of 0.01329 is calculated. As a result, the average Au concentration contributed by each nugget in these samples is 0.01329 opt, and the magnitude of *k* describes the relative size of the ‘nugget effect’ in the samples (large *k* indicates large nuggets, whereas small *k* indicates small nuggets). As a result, *k* is a fundamental sampling parameter that describes the coarseness of the gold in the samples. The value of *k* can be used in Equation 15 to define a measurement–error model for these gold results: *σ*_{tc}= , and this function is presented on a modified Thompson–Howarth graph (Stanley 1998, 2003*a*; Fig. 6).

With this fundamental sampling parameter, the associated variance stabilizing transformation that assumes only Poisson sampling error can be derived. Substituting Equation 15 into Equation 22, we obtain, after substituting the maximum likelihood value of *k*: (25) and transformation using this function produces a new distribution with stabilized (homoscedastic) variance. The original mean, standard deviation, skewness and kurtosis of these quadruplicate mean gold concentrations are 0.0524 opt, 0.0927 opt, 5.695 and 47.550 (Fig. 5), respectively, whereas the transformed mean, standard deviation, skewness and kurtosis are 3.1750 opt, 2.3890 opt, 1.741 and 5.643, respectively ( Fig. 7). As a result, this transform significantly alters the frequency distribution of the data by substantially reducing the skewness, and leads to recognition of multiple modes in the frequency distribution (Fig. 7), which may have geological or mineralogical significance. More importantly, this transformation stabilizes the sampling variance. A graph plotting the original concentration against the transformed values ( Fig. 8) illustrates the nature of this error variance transformation.

#### Binomial sampling and proportional analytical error model

Another example of the derivation and use of a variance stabilizing transformation is presented below using the 247 soil grid Cu concentrations from the Daisy Creek stratabound Cu–Ag prospect in Montana (Stanley 1984; Fig. 9). A modified Thompson–Howarth evaluation (Stanley 2003*b*) using 90 laboratory replicates with average concentrations ranging from 4 to 869 ppm, and five reference materials (three of them blind to the laboratory) with average concentrations of 4, 30, 35, 163 and 699 ppm analysed 17, 29, 21, 27 and 25 times, respectively, produce a proportional analytical error model with regression estimated *m* and *b* values of 11.3235% and 0 ppm ( Fig. 10; diamond symbols).

Calculating the total error variances of 60 additional site replicates (Fig. 10; triangle symbols) from the survey grid, and subtracting the appropriate analytical error variances, as determined by the Thompson–Howarth results for the analytical replicates (Fig. 9), yields estimates of the individual sampling error variances on these 60 site replicates. The simple linear regression model from Equation 16 (*y*=*rx*, where , *x*=−*c*^{2}+*c* and *r*=1/*n*) can then be applied to these data and optimized (*r*=Σ*x*^{2}/Σ*xy*), resulting in a value of *r*=9.48187× 10^{−6}, and a maximum likelihood value for the sampling parameter *n* of 105,464.

Note that some (25 of 60) of the calculated sampling variances for the replicate samples are negative. This is because some of the estimated analytical error variances are larger than the estimated total error variances. These negative variance values are not strictly numerically valid, but are the best estimates of the sampling error in those samples, and simply reflect the fact that in these replicates, the total error was sometimes randomly under-estimated and/or the analytical error was sometimes randomly over-estimated; note that the opposite is also likely true for some of the replicates with positive variances. As a result, all of the replicates must be used to obtain the best estimate of the sampling parameter *n* using Equation 16.

The value of *n=*105 464, along with the *m* value from the Thompson–Howarth analysis, above, can thus be used to describe a measurement–error model that involves proportional analytical error and binomial sampling error ( Fig. 11). The parameters *n* and *m* can thus be used in Equation 24 to obtain an error variance stabilizing transform for the Daisy Creek soil survey Cu concentrations. The value of *m*^{2}*n* for the Daisy Creek soil survey Cu concentrations is 1352.28, indicating that analytical error dominates in these samples (*m*^{2}*n* > 1). As a result, according to Equation 24, the appropriate error variance stabilizing transformation is: (26) because *m*=0.113235 and *n*=105 464 (Table 1, Equation T10). Applying this transformation to the soil Cu concentrations results in a distribution with skewness and kurtosis of 2.216 and 4.955, respectively ( Fig. 12). As a result, like the transformation above, this error variance stabilizing transformation produces a new variable that is not normally or near-normally distributed. Nevertheless, the transformation does produce a new variable where the propagated measurement error (as defined by the measurement–error model used in Equation 24) is constant (homoscedastic). Evaluation of this transformed data can then take place on a ‘level playing field’. A plot of the original and transformed data ( Fig. 13), analogous to Figure 8, allows comparison of this proportional/binomial error variance stabilizing transformation and the Poisson error variance stabilizing transformation.

The above error variance stabilizing transformation involves two parameters (*b* and *n*). Geochemical concentration data that are similarly subject to proportional analytical and binomial sampling errors may exhibit different values of *b* and *n*. As a result, Figures 14 and 15 present comparisons of the resulting error variance stabilizing transforms with the original concentration data for a range of parameter values. Figure 14 fixes proportional analytical error at 10%, and varies the binomial sampling parameter (*n*; Stanley 2003*a*) from 2000 to 200 000, such that different levels of binomial sampling error are accommodated by the transformation. Clearly, when sampling error is small (with large *n*), the resulting transformation produces significantly different values than when sampling error is large (with small *n*).

Similarly, Figure 15 fixes the binomial sampling parameter (*n*; Stanley 2003*a*) at 100 000, and varies the proportional analytical error from 2 to 50% such that different levels of analytical error are accommodated by the transformation. On this figure, the impacts of different amounts of analytical error are much smaller than for Figure 14, indicating that sampling error exerts far more control on the form of the resulting error variance stabilizing transformation than analytical error.

A comparison of the total error variance stabilizing transformation (Equation 26) with component analytical (logarithmic; Equation 1) and sampling (angular; Equation 3) error variance stabilizing transformations can be made by comparing Figures 16 to 19 with Figures 12 and 13. Figures 16 to 19 present the resulting histograms and transformation curves for these two component distributions, respectively. The total error variance stabilizing transformation has a form that is intermediate between the form of the component logarithmic and angular transformations. Furthermore, the angular transformation distribution most nearly approximates the distribution of the total error variance stabilizing transformation, in spite of the fact that binomial sampling error is dominant over proportional analytical error based on the value of *m*^{2}*n* and the total error variance stabilizing transformation has a logarithmic form.

Presentation of the above total error, proportional error and binomial error variance stabilizing transformation results in geographic space may lead to recognition of spatial patterns that can be associated with geological features of interest. Stanley (2006) discussed a method to present the transformed geochemical concentrations in space using a bubbleplot. As a result, Figure 20 presents a bubbleplot which depicts the original raw data, with bubble diameters proportional to the Daisy Creek soil survey Cu concentrations after rescaling to the interval 0↔1, as recommended by Stanley (2006).

Figures 21 to 23 present alternative bubbleplots depicting the total error, proportional error and binomial error variance stabilizing transformation results from the Daisy Creek soil survey Cu concentrations (with values proportional to bubble diameters). These bubbleplots depict the transformed geochemical data in different ways, and the differences between these patterns and that of Figure 20 are functions of the different transformations employed.

When comparing these bubbleplots, one should resist from founding any judgement regarding the relative success or acceptability of a transformation based on the coherence of the pattern observed. Rather, judgement should be based on the characteristics that the various transformations imparted to the transformed data (in this case, that of homoscedasticity). Nevertheless, there is a remarkable correspondence between the sizes of the bubbles in Figures 21 to 23 and the two Cu-mineralized zones on the soil grid, suggesting that these transformations do a better job of illustrating the mineralogical controls in the data than the raw data (Fig. 20).

#### Hypergeometric sampling and affine analytical error model

A last example of an error variance stabilizing transformation considers S concentrations in drillcore samples collected at a sufficient sample density relative to sample size to ensure that the binomial distribution cannot be employed to model sampling error. This example involves a lithogeochemical study of the host rocks to the Halfmile Lake South Deep Zone volcanic-hosted massive sulphide (VHMS) deposit in the Bathurst mining camp in New Brunswick (Mireku 2001; Mireku & Stanley 2005,⇓; Fig. 24). Approximately 5000 m from nine drillcores were sampled at a spacing of approximately one sample per 25 m. Samples consisted of whole NQ core, and averaged 25 cm in length (corresponding to approximately 2 kg samples). As a result, the *d*=*N*/*n* ratio for these samples is 100, a value that is low enough to ensure that the binomial distribution is inappropriate for use in estimating the sampling error.

Two pulverized in-house reference materials (the graphitic Cultus Lake shale and glacial Quadra sand, both from British Columbia) were analysed five times to monitor analytical accuracy. These, along with 20 pulverized rock samples from the survey analysed in duplicate, were also used to evaluate analytical precision. A modified Thompson–Howarth analysis of these replicate determinations (Stanley 2003*b*) reveals that S concentrations in the host rocks exhibit an affine measurement–error relationship with a relative error term of 1.7584% and absolute error term of 0.0076 wt% ( Fig. 25).

In addition to these analytical replicates, 15 sets of duplicate samples and one set of triplicate samples were collected as site replicates within 10 m of each other in the same lithology to assess local sampling error. The appropriate hypergeometric regression model is virtually identical to that of the proportional/binomial case history, above: *y*=*rx*, where , *x*=−*c*^{2}+*c*, but *r*=(*d*−1)/(*nd*−1) and *d*=*N*/*n* (Equation 27).

The resulting maximum likelihood *r* value equals 0.008276, and because *d* equals 100, the sampling parameter *n*=119.64. Thus, the appropriate hypergeometric sampling and affine analytical error model (Fig. 26): (27) defines an error variance stabilizing transformation of: (28) which after evaluating the integral becomes: (29)

With *m*=0.017584, *b*=0.0076, *n*=119.64 and *d*=100, [*m*^{2}−(*d*−1)/(*nd*−1)]=−0.0079663, so the appropriate error variance stabilizing transformation formula for the S concentrations from the wallrocks to the Halfmile Lake VHMS deposit is: (30)

A histogram of the resulting error variance stabilizing transform values reveals only a marginal decrease in skewness and kurtosis, relative to the original data (2.079 and 3.446 versus 2.538 and 5.943, respectively) indicating that the transformation did not convert the data into a near-normal distribution ( Fig. 27). This transformation, like the others above, is distinctly non-linear and concave downwards, because although the measurement–error model (Equation 30) increases with concentration ( Fig. 28), its slope shallows with concentration, indicating that the errors are naturally stabilizing at high concentrations.

Figures 29 and 30 present bubbleplots of a cross-section of the raw S concentrations and the error variance stabilizing transformation (Equation 30) through the Halfmile Lake VHMS deposit, Bathurst mining camp, New Brunswick (Mireku 2001; Mireku & Stanley 2005,⇓). These illustrate the difference between the raw and transformed results, and demonstrate that the lateral continuity of S anomalies along strike from massive sulphide mineralization is stronger in the transformed data (Fig. 30) because the intermediate S concentrations (2.5 to 5 wt% S) are emphasized. Nevertheless, the use of this error variance stabilizing transformation is not justified by its interpretability, but rather by the homoscedasticity that it creates in the transformed data.

## Discussion

Error variance stabilizing transformations for a variety of alternative hybrid models for combined sampling and analytical error can be developed using the above strategy. The set of possible error model combinations considering fixed, proportional and affine analytical error models and Poisson, binomial and hypergeometric sampling error models are presented in Table 1. Using quality control replicate analysis data, the required parameters necessary to define the appropriate alternative transformations can be estimated to allow use of a transformation that produces a new variable in which the effect of measurement error is constant. This can significantly improve a geochemist's ability to interpret geochemical data, and provides the geochemist with a specific motivation to tractably examine geochemical data from another perspective.

Other examples of error variance stabilizing transformations applied to regional geochemical survey data are presented in Stanley *et al*. (2003); Moore *et al*. (2003),⇓. Many of these make no attempt to quantify the magnitude of measurement error using replicate quality control data. Rather, these examples assume error is contributed by sampling, analysis of both, and employ the appropriate simple, unparameterized transformations (e.g. Equations 1 to 3) to extract more information from the data. This approach represents an acceptable short-cut, and can be employed routinely during geochemical data evaluation when sufficient quality control information is lacking.

The reader should also note that if the sample preparation error estimates are available and important, analogous and more complex error variance stabilizing transformations involving sampling, preparation and analytical errors (summed via their variance functions) may be developed using the above procedure.

## Conclusions

Transformation of concentration data into a new variable which exhibits a normal, near-normal or even symmetric frequency distribution is not generally possible because of the multi-modal character of many geochemical datasets. As a result, transformation for this purpose is generally misguided and unlikely to be possible.

However, transformation of concentration data to obtain new variables that exhibit other optimal properties is warranted on a variety of theoretical grounds. Use of a power transform to maximize the variance of concentration data that has been rescaled to the interval 0 ↔ 1 is: (1) theoretically achievable, (2) can be undertaken using a simple spreadsheet program, and (3) provides a resulting variable that exhibits the maximum geochemical contrast, and thus affords the geochemist with the best opportunity to recognize previously unrecognized information in the data.

Similarly, use of an error variance stabilizing transform on concentration data to create a new variable in which the effect of measurement error is constant across the range of data is: (1) theoretically achievable, (2) can be undertaken when a modest amount of quality control replicate data is available to quantify the magnitudes of different types of measurement error in the dataset, and (3) provides a result that is not biased by variable amounts of measurement error across the range of concentrations observed, thus simplifying subsequent data analysis and providing the geochemist the opportunity to recognize previously unrecognized information in the data.

In both the maximum data variance and error variance stabilizing transforms, the nature of the transformed variable's resulting frequency distribution is not a criterion used to assess whether the transform is appropriate or not. This is because the desired characteristic sought, that of maximum variance or homoscedastic measurement error, is achieved *a priori* through use of the associated transform. This affords geologists with an opportunity to avoid the temptation of employing circular logic to justify their data analysis. The data analysis procedure employed (a transformation) is not based on the interpretability of the results of the data analysis procedures, and thus is scientific. As a result, use of a maximum variance or error variance stabilizing transform will result in a more objective evaluation of geochemical data, and provide new opportunities to recognize new characteristics in the data that have not been previously been identified. Evaluation of the frequency distributions of geochemical data should routinely also consider examination of the results of these two data transformation procedures, and the results of these transformations may be effectively represented using bubbleplots with diameters that are proportional to the transformed values.

## Acknowledgements

This paper benefited from helpful discussions with Dr John Petkau, Department of Statistics, University of British Columbia, and consultations with Renjun Ma and Shou-Jun Luo of the Statistical Consulting and Research Laboratory (SCARL), Department of Statistics, University of British Columbia, as many as 14 years ago, and more recently with Drs Richard Karsten, and Jeff Hooper, Department of Mathematics and Statistics, Acadia University. It also benefited from helpful manuscript reviews by Dr A.J. Sinclair and an anonymous reviewer. It has been supported by an NSERC Discovery Grant to the author.

- © 2006 AAG/The Geological Society of London