Chapter 12
Beyond Fisher statistics

Suggested Reading
For background:
Chapter 1: Borradaile (2003)
Efron, B. and Tibshirani (1982)
Tauxe, L. et al. (1990)
To learn more:
Chapters 2-5: Fisher, et al. (1987)

12.1 Introduction

Paleomagnetists have depended since the 1950’s on the special statistical framework developed by Fisher (1953) for the analysis of unit vector data. The power and flexibility of a variety of tools based on Fisher statistics enables quantification of parameters such as the degree of rotation of a crustal block, or whether the geomagnetic field really averages to a geocentric axial dipole independent of polarity. These tools, however, require that the paleomagnetic data belong to a particular parametric distribution - the Fisher distribution.

In many important situations, the Fisher distribution fails to represent paleomagnetic data adequately. To begin with, it is likely that the geomagnetic field itself produces directions that are far from Fisher distributed (Creer et al., 1959). Collections of paleomagnetic directions from the last five million years from around the globe (e.g., McElhinny and McFadden, 1997) suggest that data from equatorial regions (Figure 12.1a) are in fact significantly elongate while those from high latitudes are more symmetric (Figure 12.1c).

The failure of the Fisher model for directional data sets makes sense if the geomagnetic field is essentially controlled by a geocentric axial dipole with “noise” added by other terms, a topic we will cover in more detail in later lectures. Such models generate spherically symmetric distributions of the VGPs. When converted to equivalent directions are more elongate as the observation site approaches the equator (see Figure 12.1b,c). Because VGPs that are farther from the pole are associated with weaker field strengths in collections of paleomagnetic data and in many models of the field, the Fisher assumption of unit vector length over emphasizes the importance of the “outliers” and leads to mean inclinations that are shallower than the true mean (see e.g., Creer, 1983).

Another example of the inadequacy of the Fisher distribution is the fact that the magnetic field exists in two stable polarity states. Because the Fisher distribution allows only uni-modal data, bi-polar data must be separated into separate modes or one mode must be “flipped” to the antipode prior to calculating a mean. Furthermore, remanence vectors composed of several components tend to form streaked distributions. Similarly, structural complications (e.g., folding) can lead to streaked distributions of directional data.

Thus, non-Fisherian data are a fact of paleomagnetic life. The Fisher-based tests can frequently be inappropriate and could result in flawed interpretations. In Lecture 11we learned the basics of Fisher statistics and how to test data sets against a Fisher model. In this lecture, we will discuss what to do when Fisher statistics fail. We will begin with essentially parametric approaches that treat certain types of non-fisherican data. We then turn to the use of non-parametric methods such as the bootstrap and jackknife in paleomagnetic applications.

12.2 Non-fisherian approaches to paleomagnetic vectors

12.2.1 The Kent Distribution

The D',I' data from equatorial sampling sites from the last five million years shown in Figure 12.1a have a more elliptical distribution than the symmetrical distribution required for a Fisherian data set. To treat such data, it is probably inappropriate to use a Fisher confidence ellipse and a distribution that allows data with elliptical directional dispersion would be better. The elliptical analogue of the Fisher distribution is the Kent distribution (Kent, 1982) and is given by:

,

If we were to collect data from the equatorial region, we might well obtain a set of directions such as those shown in Figure 12.1b. AThe Fisher ₉₅ circle of confidence for this data set is shown as the light blue line and the Kent ellipse as the heavy red line. The Kent ellipse clearly represents the the distribution of data better than the Fisher ₉₅, being elongate in the same sense as the data themselves.

12.2.2 The Bingham distribution

The Kent distribution has the advantage that it can deal with elliptical data sets while the Fisher distribution cannot. However, many paleomagnetic data sets are also bi-modal and the Kent and Fisher distributions can only deal with data sets with a single polarity. It was precisely for treating bimodal, elliptical data that the Bingham distribution was developed (Bingham, 1974). The Bingham distribution is given by:

The principle drawback of the Bingham distribution is that because the orientation matrix uses the entire data set (normal and reverse) the two modes are assumed to be antipodal and to share the same distribution parameters. The question of whether normal and reverse data sets are antipodal and have the same dispersion is in fact a question we may wish to ask! One could separate the two modes prior to calculation of the Bingham ellipse, but then the rationale for using the Bingham distribution is lost.

12.2.3 The Bingham-Le Goff approximation

Estimating the parameters for the Bingham ellipse exactly is computationally taxing and all of the available “canned” programs use the look up table of Mardia and Zemloch (1977; see Appendix). Le Goff et al. (1992) suggested some approximations which may be valid for concentrated distributions. They also introduced the concept of weighting results according to some reliability criteria. For the general case, however, it seems preferable to use the exact Kent (1982) ellipses on uni-modal data sets. These could of course be weighted if such weighting is desired.

12.2.4 The bi-Gaussian Distribution

Until now we have continued the Fisher assumption of unit vectors. As already mentioned, neglect of the vector strength can lead to bias. Love and Constable (2003) began the hard work of incorporating intensity information into the parameter estimation problem. Their method can handle bi-modal spherical Gaussian data such as those shown in Figure 12.3. Estimation of the Love parameters are beyond the scope of these lecture notes. Moreover, many data sets are not spherically symmetric as already noted and the Love and Constable (2003) approach must be generalized to elliptical, more “blade-like” data sets than the “cotton balls” currently treatable.

Figure 12.4: Hypothetical non-Fisherian data set. Normal and reversed polarity data that are not symmetrically distributed. Filled (open) circles plot on the lower (upper) hemisphere. b) Equal area projection of 500 bootstrapped means for para-data sets drawn from the data shown in a) Tiny circles (plusses) are lower (upper) hemisphere projections from the normal (reversed) polarity mode of each para-data set. c) Normal polarity data from a) transformed into VGPs and plotted in polar equal area projection. The dashed ellipse is obtained from VGP transformation of 95, while the solid ellipse is determined by a bootstrap.

12.3 The simple bootstrap

As we have mentioned, “real” data may be pathological in several respects including bi-modal and elliptically distributed data. Neither the Kent nor the Bingham statistical methods have the test for common mean so critical to paleomagnetic studies. Moreover, none of the methods we have described so far can provide confidence ellipses for an off-center mean direction as is likely to occur in records of the geomagnetic field (see Figure 12.1c). Finally, data may be overprinted or contain the record of a paleomagnetic transition, resulting in “streaked” or non-antipodal mean directions, conditions that make the conventional methods inappropriate. In this section we will discuss alternative methods for estimating confidence bounds which is sufficiently flexible to accomodate all of these short comings, provided the data set is large enough.

In Figure 12.4a we show a fairly typical “not great” paleomagnetic data set. The data are elliptical, bi-modal and one has the suspicion that the normal and reverse modes may be neither antipodal nor share the same concentration or ovalness parameters. Clearly some non-parametric approach would be desirable. The approach for characterizing uncertainties for vectors we will take here is based on a technique known as the statistical bootstrap. As we shall see, the bootstrap has the flexibility to allow us to treat awkward data sets like that shown in Figure 12.4a.

The principles of the bootstrap are giving in the Appendix. In essence, the parameter of interest (say, the mean vector) is calculated from many resampled data sets, whose data points are selected at random from the original data. The bootstrapped estimates “map out” the likely distribution of the parameter, allowing estimation of confidence regions. Before we extend the bootstrap from the scalar treatment in the appendix to vectors, it is important to point out that with the bootstrap, it is assumed that the underlying distribution is represented by the data, demanding that the data sets be rather large. Moreover, the bootstrap estimates are only asymptotically true, meaning that a large number of bootstrap calculations are required for the confidence intervals to be valid. It’s a good thing we have fast computers with huge hard-drives.

There are a variety of ways we can use the bootstrap to estimate confidence regions for paleomagnetic data. We will start with the most “Fisher” like approach of taking unit vectors of a single polarity. Then we will accommodate dual polarity data sets and develop analogous tests to those so useful for Fisher distributions.

To do a simple bootstrap on a data set with only one polarity (say the normal data in Figure 12.4a, we first randomly draw N data points from the data shown in Figure 12.4a. Each set of N data points is a “para-data set” (also known as a “pseudo-sample”). We then calculate a Fisher mean of the para-data set (one little circle in Figure 12.4b. This resampling procedure can be repeated many times. We show 500 such bootstrapped means in Figure 12.4b.

Now we can estimate the region of 95% confidence for the bootstrapped means. A non-parametric approach would be to draw a contour enclosing 95% of the bootstrapped means. In many applications, paleomagnetists desire a more compact way of expressing confidence regions (for example, to list them in a table) and this necessitates some parametric assumptions about the distribution of the means. For this limited purpose, approximate 95% confidence regions can be estimated by noting the elliptical distribution of the bootstrapped means and by assuming that they are Kent (1982) distributed.

When paleomagnetic data are bimodal, as in Figure 12.4a, we can proceed in one of two ways. We could just calculate the principle eigenvector of the orientation matrix (V₁) as in Bingham statistics of each bootstrapped para-data set or we can separate the data into two modes and calculate Fisher means for each mode separately (as in Figure 12.4b).

To separate the data into normal and reverse sub-sets, we first calculate the orientation matrix and its principal eigenvector V₁. By transforming the data set into the coordinates of V₁ (see Appendix to Lecture 2), the transformed inclinations could be used to assign polarity in an automated and unbiased way; positive inclinations belong to one mode and negative inclinations to the other. After separation, Fisher means can be calculated for each mode separately. Alternatively, if a more robust estimate of the “average” direction is desired, one could calculate the principal eigenvector V₁ of each mode, which is less sensitive to the presence of outliers.

Figure 12.5: The bootstrap test for a common mean. a) Equal-area projections of two data sets two simulated Fisherian data sets (triangles and circles) each with of 20. b) Means and 95s of data sets shown in a). c) Histograms of x1 components of the bootstrapped means from 500 para-data sets. Also shown are the bounds for each data set that include 95% of the components. The confidence intervals for the different data sets overlap. d) Same as c) but for the x2 component. The confidence intervals do not overlap. e) same as d) but for the x3 component.

12.4 The parametric bootstrap

The bootstrap just described is a “simple” or “naive” bootstrap in that no distribution is assumed for the data. We did assume, however, that all the uncertainty inherent in the data is reflected in the data distribution. If the data set is smaller than about N = 20, this leads to uncertainty ellipses that are too small (Tauxe et al., 1991). Many paleomagnetic data sets are smaller than this, yet they are demonstrably non-Fisherian. Luckily, if we are able to assume some parametric form for data from e.g., a given site, we can perform a superior technique which is known as the parametric bootstrap. As applied here, we assume that each site with N_s samples is Fisher distributed (in principle, a testable assumption). Then, after random selection of a particular site for inclusion in the para-data set, we draw N_s new directions from a Fisher distribution with the same mean direction, and N. From these, we then calculate a substitute mean direction, and use that in the para-data set. Otherwise, we follow the same procedure as before.

For large data sets (N > 25), the parametric and simple bootstraps yield very similar confidence ellipses. For smaller data sets, the parametric ellipses are larger, and are probably more realistic.

12.5 Application to VGPs

We stated earlier that the transformation of Fisher distributed directional data to VGPs can cause a distortion from rotationally symmetric data to elliptically distributed data. To illustrate what happens if data are distributed as in our hypothetical normal data set, we calculate VGPs in Figure 12.4a, the ellipse described by dp,dm (Lecture 11), and our bootstrap confidence ellipses. The VGPs are plotted in Figure 12.4c and the dp,dm is shown as a dotted line. Because the declinations are smeared and the dp must point towards the site (shown as a triangle), the long axis of the so-called 95% confidence regions is perpendicular to the actual data distribution. The 95% confidence ellipse, calculated using the bootstrap, is shown as a solid line. The bootstrapped confidence ellipse gives a better sense of the uncertainty and follows the trend in the data.

12.6 When are two data sets distinct?

The test for a common mean addresses the question “can the means of two data sets be discriminated from one another?” Another way of putting it is, “If a set of bootstrap means is examined, are there two distinct groups or is there just one?” We explore these ideas by considering the same Fisherian data sets we used in Lecture 11for the Watson’s V _w test. In Figure 12.5 we show two data sets (triangles and circles), each drawn from distributions with a of 20. The mean direction of each lies outside the confidence region of the other and the but the V _w test of Watson has a value of 11.7 with a critical value of 6.3, hence they fail the test for a common mean.

In order to develop a bootstrap test analagous to the V _w test for use on non-Fisherian data sets, we first convert a set of bootstrapped mean directions to to cartesian coordinates. Histograms of the cartesian coordinates of the bootstrap means (Figure 12.5c-e) are distinct in the x₂ component, confirming that the two means can be distinguished at the 95% confidence level.

12.7 Application to the “reversals test”

The so-called reversals test in paleomagnetism constitutes a test for a common mean for two modes, one of which has been “flipped” to its antipode. We apply our bootstrap test for common mean to the data shown in Figure 12.4. The histograms of the cartesian coordinates of the bootstrapped means are shown in Figure 12.6. There are two “humps” in the bootstrap test. However, the confidence intervals for the normal and reverse antipodes overlap, thereby suggesting that the two means cannot be distinguished at the 95% level of confidence. Thus, the data in Figure 12.4 pass the bootstrap reversals test.

12.8 Application to the fold test

A final test is extremely useful in paleomagnetism: the fold test (Lecture 9). One of the key components in paleomagnetic studies is to determine the coordinate system (geographic, tilt adjusted or somewhere in between) for which the directional data are most tightly clustered. If a rock has moved from its original position, was it magnetized in the original, in the present or in some other position? Moreover, is simple rotation about strike an appropriate method to restore the beds to their original positions? In the classic fold test as first proposed by Graham (1949), the directions of magnetization of a deformed rock unit are assumed to be most closely parallel in the orientation in which the magnetization was acquired. Therefore, if a rock has retained an original magnetization through a subsequent folding or tilting event, the magnetic directions may cluster most tightly after they have been rotated back to their original positions. This of course is NOT true for elongate data such as those shown in Figure 12.1a.

Figure 12.8: The largest eigenvalues 1 of the orientation matrices from representative para-data sets drawn from Figure 12.6. The directions are adjusted for tilt incrementally from -50% to 150%. The largest value of 1 occurs near 100% in all of the para-data sets. Histogram is of 500 maxima of 1 and the calculated 95% confidence interval. These data “pass” the bootstrap fold test. [Redrawn from Tauxe and Watson, 1994.]

The fold test appears at first glance to be simple, but it is not. The primary problem is that paleomagnetic vectors are never perfectly parallel. The scattered nature of the data means that a statistical test is necessary to determine whether clustering is “significantly” better in one orientation or another.

In Lecture 11we suggested that variances could be compared using an F-test, so it was long the practice in paleomagnetism to compare estimated precisions before and after tilt adjustment (McElhinny, 1964). The ratio of the two estimates of were compared with those listed in statistical “F” tables (see Appendix to Lecture 11. Ratios higher than the “F” value for a given N were deemed to constitute a significant increase in concentration after adjusting for tilt, thus representing a positive fold test. This test can be done on the back of an envelope and is still in frequent use.

Although its simplicity is a great strength, there are several problems with the classical fold test. First, the geomagnetic field has two preferred states and is not perfectly dipolar. Directions observed in paleomagnetic samples are therefore not only scattered but are often of two polarities. Second, the magnetic directions may be most tightly clustered somewhere other than in “geographic” or 100% tilt adjusted coordinates. Finally, structual “corrections” are not perfectly known. Not only are the bedding orientations themselves often difficult to measure accurately, but detection of complications such as plunging folds, and multiple phases of tilting requires extensive field work. It is nearly impossible to assess rotation about the vertical axis on the basis of field relations alone, as it results in no visible effect on the dip of the beds themselves. Because of this uncertainty, we might reasonably ask whether if the data are actually most tightly clustered at, say 90% tilt adjusted (as opposed to 100%), does this constitute a “failed” fold test.

We consider first the problem of dual polarity. We plot a data set in geographic coordinates in Figure 12.7a and in tilt adjusted coordinates in Figure 12.7b. The polarity is ambiguous and the calculation of necessitates using the tilt adjusted data to identify the polarity of the samples. The classic fold test requires calculation of which can only be done with data of a single polarity. Obviously, fold tests that rely on will not be straightforward with data such as these.

An alternative approach is based on the orientation matrix. In the orientation matrix, polarity does not play a role and the “tightness” of grouping is reflected in the relative magnitudes of the eigenvalues (). As the data become more tightly grouped, the variance along the principal axis grows and those along the other axes shrink. Thus, examination of the behavior of ₁ during unfolding would reveal the point at which the tightest grouping is achieved, without knowledge of polarity.

Suppose we find that the degree of unfolding required to produce the maximum in ₁ is 98%. Is this a positive fold test suggesting a pre-folding remanence or is the difference between 98% and 100% significant? For this we call on the familiar bootstrap. Numerous para-data sets can be drawn. We can then calculate the eigenparameters of the orientation matrix for a range of percent unfolding. Some examples of the behavior of ₁ during tilt adjustment of representative para-data sets drawn from the data in Figure 12.7a are shown in Figure 12.8. Figure 12.8 is a histogram of maxima of ₁ from 500 para-data sets. These are sorted and the 95% confidence interval for the degree of unfolding required to produce the tightest grouping (the highest ₁) is thus constrained to lie between 97 and 102%.

The data from Figure  12.7a are shown after 100% tilt adjustment in Figure 12.7b. The tilt adjusted data are not only better grouped, but now the polarities of most samples can be readily determined. An advantage of the bootstrap approach is the fact that the data do not need prior editing to split them into normal and reversed polarity groups, which is a particularly onerous task for the data considered here.

For small data sets, we could of course employ a parametric bootstrap, whereby para-data sets are generated by first randomly selecting a site for inclusion, then by drawing a substitute direction from a Fisher distribution having the same D, I, N, and .

Appendix

A Estimation of Kent 95% confidence ellipse

Kent parameters are calculated by rotating unimodal directions x into the data coordinates x' by the transformation:

(A1)

where = (₁,₂,₃), and the columns of are called the constrained eigenvectors of T. The vector ₁ is parallel to the Fisher mean of the data, whereas ₂ and ₃ (the major and minor axes) diagonalize T as much as possible subject to being constrained by ₁ (see Kent (1982), but note that his x₁ corresponds to x₃ in conventional paleomagnetic notation). The following parameters may then be computed

(A2)

As defined here, = R/N (R is closely approximated by the equation for R in Lecture 11). Also to good approximation, ₂² = ₂, and ₃² = ₃, where _i are the eigenvalues of the orientation matrix. The semi-angles ₉₅ and ₉₅ subtended by the major and minor axes of the 95% confidence ellipse are given by:

(A3)

where g = -2 ln(0.05)/(N²).

The tensor is, to a good approximation, equivalent to V, the eigenvectors of the orientation matrix. Therefore, the eigenvectors of the orientation matrix V give a good estimate for the directions of the semi-angles by:

(A4)

where for example the x₂ component of the smallest eigenvector V₃ is denoted v₂₃.

B Estimation of Bingham 95% confidence parameters

The Bingham distribution is given by:

To estimate the axes of the Bingham confidence ellipse, we first calculate the eigenparameters of the orientation matrix as for Kent parameters and described in the Appendix to Lecture 9. Note that through-out these lecture notes _i are the eigenvalues and V_i are the eigenvectors. The principle eigenvector V₁ of the orientation matrix is associated with the largest eigenvalue ₁. In Bingham (1974), ₁ is the ₃ and ₃ is ₁. In Bingham statistics, the V₁ direction is taken as the mean - it is not always parallel to the Fisher mean of a unimodal set of directions.

The maximum likelihood estimates of k₁,k₂, the concentration parameters are gotten by first maximizing the log likelihood function:

Bingham (1974) set k₃ = 0, so the semi axes of the confidence ellipse about the principle direction V₁, associated with ₃, are therefore:

Because k₁ < k₂ < 0, the semi-axes are positive numbers. Please note that we use the corrected version of Tanaka (1999) as opposed to the more oft quoted but erroneous treatment of Onstott (1980) in this appendix. Note also that the N is required for because we have normalized the s to sum to unity for consistency with other eigenvalue problems in this set of lecture notes. The N is missing in the treatment of Tanaka (1999) presumably because the eigenvalues sum to N. Finally, note that these values of are in radians and must be converted to degrees for most applications.

C The statistical bootstrap

In Figure B1, we illustrate the essentials of the statistical bootstrap. We will develop the technique using data drawn from a normal distribution. First, we generate a synthetic data set by drawing 500 data points from a normal distribution with a mean of 10 and a standard deviation of 2. The synthetic data are plotted as a histogram in Figure B1a. In Figure B1b we plot the data as a Q-Q plot against the z_i expected for a normal distribution (see Abramowitz and Stegun (1970).

In order to calculate the appropriate values for z_i assuming a normal distribution:

1. For i = 1 N, calculate p = i __ N+1.
2. If p > 0.5, then q = 1 - p; if p < 0.5, then q = p.
3. Calculate the following for all p0.5:
   

   and

   where a₁ = 2.515517,a₂ = 0.802853,a₃ = 0.010328,a₄ = 1.432788,a₅ = 0.189269,a₆ = 0.001388.
4. If p > 0.5, then z_i = u; if p < 0.5, then p = -u.
5. If p = 0.5, then z_i = 0.

The values of z_i calculated in this way for the simulated Gaussian distribution are plotted as the “normal quantile” data in Figure B1b.

The simulated data fall along a line in Figure B1b which suggests that they are normally distributed (as one would hope). To test this in a more quantitative way, we can calculate D_N⁺ and D_N^- as follows:

1. Calculate the mean and standard deviation for the data.
2. Then calculate:

   and

   
3. Substitute q into the following expression (function 7.1.26 from Abramowitz and Stegun (1970):
   where a₁ = 0.254829592,a₂ = -0.284496736,a₃ = 1.421413741, and a₅ = 1.061405429.
4. Change the sign of erf(q) such that it has the same sign as q.
5. Substitute F(x) = 0.5(1 + erf(q)) into Equations D3 and D4 in the Appendix to Lecture 11for D_N⁺ and D_N^- respectively. The Kolmogorov-Smirnov parameter D (e.g., Fisher et al., 1987) is the larger of D_N⁺ or D_N^-.
6. The null hypothesis that a given data set is normally distributed can be rejected at the 95% level of confidence if D exceeds a critical value D_c given by 0.886/.

Applying the foregoing to the data in Figure B1a yields a D value of 0.0306. Because N = 500, the critical value of D, D_c at the 95% confidence level is 0.0396. Happily, our program has produced a set of 500 numbers for which the null hypothesis of a normal distribution has not been rejected. The mean is about 10 and the standard deviation is 1.9. The usual Gaussian statistics allow us to estimate a 95% confidence interval for the mean as ±1.96/ or ±0.17.

The mean of the synthetic dataset is about 10 and the standard deviation is 1.9. The usual Gaussian statistics allow us to estimate a 95% confidence interval for the mean as ±1.96/ or ±0.17.

In order to estimate a confidence interval for the mean using the bootstrap, we first randomly draw a list of N data by selecting data points from the original data set. Some data points will be used more than once and others will not be used at all. We then calculate the mean of the “para-data set”. We repeat the procedure of drawing para-data sets and calculating the mean many times (say 10,000 times). A histogram of the “bootstrapped” means is plotted in Figure B1c. If these are sorted such that the first mean is the lowest and the last mean is the highest, the 95% of the means are between the 250^th and the 9,750^th mean. These therefore are the 95% confidence bounds because we are approximately 95% confident that the true mean lies between these limits. The 95% confidence interval calculated for the data in Figure B1 by bootstrap is about ± 0.16 which is nearly the same as that calculated the Gaussian way. However, the bootstrap required orders of magnitude more calculations than the Gaussian method, hence it is ill-advised to perform a bootstrap calculation when a parametric one will do. Nonetheless, if the data are not Gaussian, the bootstrap provides a means of calculating confidence intervals when there is no quick and easy way. Furthermore, with a modern computer, the time required to calculate the bootstrap illustrated in Figure B1 was virtually imperceptible.

Bibliography

 

 

[ Home ]  [ Databases ]  [ Events ]  [ Tools ]  [ Publications ]  [ Links ]  [ Google Site Search ] [ Metadata ]  [ Who's Who ]  [ Browsers and Plugins ]  [ Database Indexes ]
This page was last updated on 04-Apr-2008 Sponsored by NSF EAR 0000998	Supported by the San Diego Supercomputer Center and the Scripps Institution of Oceanography