On point estimation of the abnormality of a Mahalanobis index.

Mahalanobis distance may be used as a measure of the disparity between an individual's profile of scores and the average profile of a population of controls. The degree to which the individual's profile is unusual can then be equated to the proportion of the population who would have a larger Mahalanobis distance than the individual. Several estimators of this proportion are examined. These include plug-in maximum likelihood estimators, medians, the posterior mean from a Bayesian probability matching prior, an estimator derived from a Taylor expansion, and two forms of polynomial approximation, one based on Bernstein polynomial and one on a quadrature method. Simulations show that some estimators, including the commonly-used plug-in maximum likelihood estimators, can have substantial bias for small or moderate sample sizes. The polynomial approximations yield estimators that have low bias, with the quadrature method marginally to be preferred over Bernstein polynomials. However, the polynomial estimators sometimes yield infeasible estimates that are outside the 0-1 range. While none of the estimators are perfectly unbiased, the median estimators match their definition; in simulations their estimates of the proportion have a median error close to zero. The standard median estimator can give unrealistically small estimates (including 0) and an adjustment is proposed that ensures estimates are always credible. This latter estimator has much to recommend it when unbiasedness is not of paramount importance, while the quadrature method is recommended when bias is the dominant issue.

Introduction

The Mahalanobis distance is frequently used in multivariate analysis as a statistical measure of distance between a vector of scores for a single case and the mean vector of the underlying population or a sample of data. It was developed by Mahalanobis (1936) as a distance measure that incorporates the correlation between different scores. See also DasGupta (1993). The Mahalanobis distance of a vector , of say variables (scores), from a population mean is defined as where is the population covariance matrix. The square of the Mahalanobis distance, , is sometimes referred to as the Mahalanobis index (Huberty and Olejnik, 2006, p. 271). If the population follows a multivariate normal distribution (MVN) and is an observation from this same distribution, then the Mahalanobis index follows a central chi-square distribution on degrees of freedom. In this paper, interest focuses on estimating , the proportion of the population that gives a more unusual Mahalanobis index than where is a specified vector, under the assumption that the population distribution is a MVN distribution. That is where . For example, might be a patient’s profile from a set of medical tests, when would be the proportion of the population with a profile that is more unusual than that of the patient.

The corresponding Mahalanobis distance in a sample, of say observations, is defined as where and are the sample mean vector and sample covariance matrix, respectively. Under the assumption that and the sample data are from the same MVN distribution, the sample Mahalanobis index is proportional to a central distribution with and degrees of freedom. See, for example, Mardia et al. (1979).

We were initially motivated by the need to estimate the abnormality of a single patient’s profile in neuropsychology. The problem arises, for example, when psychologists need to assess how a patient with some brain disorder or a head injury is different from the general population or some particular subpopulation. This assessment is usually based on the patient’s scores in a set of tests that measure different traits or abilities. The abnormality of the case’s profile of scores can then be expressed in terms of the Mahalanobis index between this profile and the mean of the normative population or normative sample. The degree of abnormality is measured by where is the case’s profile and is treated as a fixed quantity.

A Hotelling’s significance test for testing whether the case could belong to the normative population is proposed in Huizenga et al. (2007). Their test is based on the central distribution to which the Hotelling’s test statistic is proportional. Crawford et al. (2016) give a confidence interval for the probability () of getting a more extreme profile than the case. The confidence interval is based on a non-central distribution with a non-centrality parameter that is proportional to the case’s Mahalanobis index. The confidence intervals are correct, in that their coverage levels equal the nominal confidence level exactly. In contrast, the -value from the Hotelling’s test provides an obvious point estimator of , but it is biased. Indeed, the problem of finding an unbiased estimator of has not been resolved.

Here we consider a number of obvious estimators of and propose some new, less obvious estimators. The bias and mean square error of all the estimators are compared in extensive simulations. No estimator is uniformly better than all alternatives, but a small selection of the estimators is clearly to be preferred. As well as bias and mean square error, other criteria and desirable qualities in an estimator are also considered. In this paper, no distributional assumptions are made about the source of , other than when testing whether could be the profile of a member of the normative population.

The need to estimate the value of for Mahalanobis distances does not only arise in psychology. In the literature, the commonly used estimates of are the -value computed from the chi-square distribution of the sample Mahalanobis index, or the -value from the central distribution associated with Hotelling’s test. For example, in remote sensing image analysis, Foody (2006) was interested in measuring the closeness of an image pixel to a single class centroid. For that, he used the Mahalanobis distance and converted the calculated Mahalanobis distance, of a particular image pixel from a specified class centroid, to its associated -value from the chi-square distribution. He then interpreted the -value as the probability of obtaining a Mahalanobis distance as extreme as that observed for a particular pixel with respect to a specified class, thus effectively equating the -value to .

In environmental and health science, Liu and Weng (2012) used Mahalanobis distance in public health studies to enhance the resolution of satellite imagery. They conducted a spatial–temporal analysis of West Nile Virus outbreak in Los Angeles in 2007 using sensing variables and infective mosquito surveillance records. Mahalanobis distance was used to identify and map the risk areas where habitat was suitable for infective mosquitoes. Liu and Weng (2012) calculated the distance between a vector of environmental variables and the mean vector of environmental factors at the closest locations of mosquito infections. Locations with smaller values of Mahalanobis distances indicated a more favorable habitat for the mosquitoes and hence an area of higher risk. They assumed that Mahalanobis distance follows a chi-square distribution, from which was calculated for each map pixel. Pixels with between 0.6 and 0.9 (0.9 and 1.0) were considered moderate risk (high risk) areas and then a risk map was produced.

In analytical chemistry, Shah and Gemperline (1990) were interested in analyzing the near-infrared reflectance spectra of raw materials. They used Mahalanobis distance as a classification technique for pattern recognition to classify new samples by comparing them to measurements of predetermined classes. Each sample was classified according to the -value associated with its Mahalanobis distance from the class centroids. Shah and Gemperline (1990) needed to estimate the -value for each new sample and used the chi-square distribution to estimate these probabilities. They considered samples with probability levels between (0–0.01), (0.01–0.05) or (0.05–1.0) to be nonmembers, outliers or members, respectively.

A sample Mahalanobis distance has an exact distribution. Hence, unsurprisingly, the distribution has also been frequently used to quantify the rarity/commonness of a Mahalanobis distance. For example, Lu et al. (2005) used the two groups Hotelling’s test for detecting differential expressions in genetic microarrays. They conducted a microarray experiment in which samples from a disease group and from a normal group were obtained. They based their test for gene differential expression on the scaled distribution of the Hotelling’s statistic. Some other important applications of the Hotelling’s statistic, Mahalanobis distance and the associated -values include, for example, multivariate outlier detection (e.g. Garrett, 1989; Hardin and Rocke, 2005) and multivariate quality control charts (e.g. Sullivan and Woodall, 1998; Johnson and Wichern, 2007). However, some methodological researches in causal inference argue that Mahalanobis distances do not work very well when the dimension, , of is greater than 8 or the normality assumption is not fulfilled. See, for example, Stuart (2010) and the references therein.

We conducted a simulation study to test the performance of the sample -value associated with the distribution, denoted by , and that associated with the chi-square distribution, denoted by , in estimating the probability . Simulation results show that both are biased estimates of . We propose some alternative estimators of and compare them in terms of their bias and root mean square error in the simulation study. Some of the proposed estimates have much lower biases than the estimators derived from the and chi-square distributions.

Three of the alternative point estimators of are based on its confidence intervals. The first uses the frequentist median of the non-centrality parameter, and is denoted by . The second proposed estimator uses the Bayesian median of the non-centrality parameter or its frequentist median, whichever is greater. We call it the modified median estimator and denote it by ; and only differ when approaches 100%. The third estimator in this group is a Bayesian estimator; it is based on the idea of probability matching priors and is denoted by . We propose another two new estimators of based on the mean of the non-centrality parameter of a non-central distribution; these are denoted by and . Estimators derived from a Taylor expansion ()–Bernstein polynomials of degree 4, 7 and 10 (, and ) and a quadrature polynomial approximation of degree 4, 7 and 10 (, and )–are also proposed and shown to be approximately unbiased in the broad range of situations examined in the simulation study.

The paper is organized as follows. In Section  2, we briefly discuss the two frequently used point estimators of . The new twelve proposed point estimators of are detailed in Section  3. All the fourteen point estimators are then compared in Section  4, where we present and discuss the results of the simulation study. In Section  5, we examine the behavior of each estimator at different observed values of the Mahalanobis index. We also briefly consider the median error of the estimators (rather than average error) and mean absolute error. Concluding comments are given in Section  6.

Two plug-in maximum likelihood estimators of

We aim to derive an unbiased estimate of where is the population Mahalanobis index , which follows a chi-square distribution on degrees of freedom, and is the Mahalanobis index of the case. That is equals , where is the vector of a case’s profile of scores from tests. The probability is the proportion of the population that has a profile more extreme than the case. A minimum variance unbiased estimator of is readily available (see Section  3.2) but obtaining an unbiased estimator of is much harder.

Let and denote the maximum likelihood estimates of and , respectively, hence . Simple estimates of can be obtained by replacing the unknown parameters in Eq. (2) with their maximum likelihood estimates. This gives our first estimator, It is well-known that where is Hotelling’s statistic and is a central distribution on and degrees of freedom. Then is the -value from testing the null hypothesis that the case is a member of the control population. Consequently, it is commonly used as a point estimate of the proportion of the normative population with more extreme profiles than the case.

Another frequently used plug-in maximum likelihood estimate, denoted by , is the -value from the chi-square distribution. Instead of replacing and by and everywhere in (2), is obtained by only making this replacement on the right-hand side of the inequality. Thus and we have that Simulation results in Section  4 show that is generally better than as an estimate of . However, both are biased and underestimates in most cases, with absolute bias that is getting higher for larger values of the true parameter .

New point estimators of

Estimators derived from confidence intervals

Classical estimator of the median

Based on the work of Reiser (2001), Crawford et al. (2016) proposed a method for constructing confidence intervals on . The observable sample statistic has a non-central distribution with , degrees of freedom, respectively, and a non-centrality parameter , i.e.  To construct a confidence interval for , define as the value of for which is the -quantile of . Then, a confidence interval for is given by where is the cdf of a chi-square distribution on degrees of freedom.

Using the same technique, a point estimator of is given by its median estimate. We have that is the value of at which is the median of the distribution. The first of our new estimators, , is defined as Although is a biased estimator of , simulation results in Section  4 show that it usually has a smaller bias and mean square error than at all values of the true parameter .

Modified estimator of the median

As decreases, so does the median of . Since , a lower bound on the median of is the median of . ( is the ordinary central distribution on and degrees of freedom.) If is less than this lower bound, one approach is to set to zero. This is the standard approach adopted in the construction of confidence intervals, where the same problem arises, as discussed in Reiser (2001). The problem arises whenever is small, even if it is above the lower bound. To illustrate, suppose , and , so . Then calculation gives . Thus when a patient’s estimated Mahalanobis index is 0.4, then 0.51% is the estimate of the proportion of the normal population with a smaller true Mahalanobis index than the case. However, if 0.4 were the true Mahalanobis distance of the case, then the actual proportion of the normal population with a smaller Mahalanobis distance than the case is calculated at 1.75% (), from a chi-square distribution on 4 degrees of freedom. The disparity between 0.51% and 1.75% is substantial and, moreover, intuitively one would expect uncertainty to result in being less extreme than , rather than being greater than it. As decreases the situation worsens. When , , while when .

A pragmatic solution was proposed by Garthwaite et al. (2016). They supposed an individual’s sample Mahalanobis index was and considered the question: “What proportion of the population will have a true Mahalanobis index that is bigger than this individual?” under two situations They argue that the answer in situation (i) should be no bigger than in situation (ii). They suggest that the proportion should be estimated for each situation and the smaller estimate selected. Adopting that approach, we construct another estimator, , as follows.

the individual is the case

the individual is a randomly chosen member of the population.

Let be the Mahalanobis index of a randomly selected individual from the population and let be the proportion of the population with a larger Mahalanobis index than . Then is a random variable and, before observing , has a uniform distribution over the interval (0, 1). In consequence, has a chi-square distribution on degrees of freedom. This chi-square distribution can be taken as the prior distribution in a Bayesian analysis in which there is a single datum, . Note that there is nothing arbitrary about this prior distribution; it is the distribution of because . The likelihood follows from Eq. (10) as . We obtain the posterior distribution of , and compute its normalizing constant through numerical integration. A simple search procedure is used to find the posterior median of , say . We then use to enhance the median estimator in (12) and propose the modified median estimator, , as

Obviously, and only differ when is small and then the differences are slight in absolute terms , though is far from 1 for very small . This is illustrated, for and , in Fig. 1, Fig. 2, where and are plotted against at different observed values of and , respectively. As takes the lower value of and , Fig. 1 shows that and are identical for , while Fig. 2 shows that, as gets small, is clearly less than 100% (as common sense dictates it should be) while approaches 100%.

Fig. 1

and at .

Fig. 2

and at .

As expected, simulation results in Section  4 show that the bias and mean square error of the modified estimator are nearly identical to those of the median estimator . We recommend over for use in practice to avoid the problem of getting , as discussed above.

Bayesian probability matching

Bayesian credible intervals quite often have the same endpoints as frequentist confidence intervals if the Bayesian intervals are based on uninformative prior distributions. Indeed, there has been substantial interest in probability matching priors  (Datta and Mukerjee, 2004), which are designed to give credible intervals that match confidence intervals. To construct our next estimate of , we suppose a prior distribution has been found that gives posterior credible intervals which match the confidence intervals specified in Eq. (11). Treating the confidence intervals as exact credible intervals, they determine a posterior distribution for the proportion . We use a sufficiently large number, say , of one-sided credible interval limits to construct the posterior distribution. The posterior mean is then used as a point estimate, say , of .

Specifically, we estimate the interval limit as the value of for which is the -quantile of , for . As in (11), the -quantile of the posterior distribution of is given by . The posterior mean is then computed as In practice, we take , so that we use the quantiles which is a sufficiently fine partition for our purpose and is not expensive in computing time.

Based on simulation results in Section  4, the estimate is a badly biased estimate of . This result illustrates an important fact: while posterior distributions obtained from exact probability matching priors will (by design) give interval estimates with good frequentist properties, the posterior mean may be far from meeting the frequentist definition of unbiasedness.

Estimators based on the mean of

Our next proposed estimators of are based on the estimated mean value of , say . If is given by Eq. (10), then is the uniformly minimum variance unbiased estimator of the non-centrality parameter of the non-central distribution . However, it is well-known that this is not always positive and is therefore inadmissible. See, for example, Johnson et al. (1995). To avoid a negative estimate of , put Using , we propose the estimator of as Unfortunately, based on the simulation results in Section  4, can have marked bias as an estimator of .

The estimator in (16) is also inadmissible (Chow, 1987), but Rukhin (1993) showed that, for , is an admissible estimator of . We base our next estimate, , of on and put However, as illustrated in Section  4, is generally better than in terms of bias and mean square error.

An estimator based on a Taylor expansion

We expand the cdf of the chi-square distribution about as where is the pdf of a chi-square distribution with degrees of freedom. We set equal to in Eq. (15) and take the expected value of both sides of (20). This gives where, from the variance of with , is given by with , and . (As defined earlier, .)

In the light of (21), to obtain an approximately unbiased estimate of to the second order, it seems natural to base an estimate on say, such that We start with the case where is greater than the mode of the chi-square distribution. In this case is negative, and we can write for some . From (23), (24) we have We define another estimate, , by replacing in (24) with : Suppose is large relative to . If , then and will be better than as an estimate of . If , then and will again be better than as an estimate of . On the other hand, supposing that is small relative to , then and . The consequence is that defined in (26) is expected to be better than , in terms of the mean square error, as an estimate of . The other case in which is less than or equal to the mode of the chi-square distribution can be treated similarly.

It remains now to estimate the right hand side of (26). We find an unbiased estimate, say , of expressed as where , and are chosen such that . Specifically, equating the corresponding coefficients of in to those in (22), we get , and . It is straightforward to show that However, no simple unbiased estimate can be found for , instead, we estimate it as . The estimator is finally expressed as Using this Taylor based estimate, our proposed approximately unbiased estimate of is given by Simulation results show that is usually one of the better estimates of . More information is given in Section  4.

Estimators based on polynomial approximations

None of the point estimators we have proposed so far attain approximate unbiasedness uniformly for all values of . This motivates another set of point estimators that are approximately unbiased uniformly for all . Using polynomial approximations, we aim to base the proposed estimator of in this section on a good global estimate of the non-centrality parameter . This means that in searching for an approximately unbiased estimate of , we cover wide areas of the chi-square cdf and do not locally search around some estimate of , as was proposed in Section  3.3 when using the Taylor expansion. In principle, estimates based on global approximation of the cdf should prove better, in terms of bias, than an estimate based on a local approximation.

We introduce a set of unbiased estimates of , denoted for now by , which are based on approximating the probability in (5) as a polynomial function of degree in . From Weierstrass’s Theorem, any function of a variable, say, can be approximated by a polynomial of , provided the function satisfies weak regularity conditions. Now is a function of that meets these regularity conditions, so we may put The coefficients are known functions in (see below).

The key to exploiting Eq. (31) is that the moments of are also polynomials in . Specifically, with defined by Eq. (10), the th moment, , is a polynomial of of degree . Writing in the polynomial form (31), it can therefore be estimated by another polynomial in as follows where the coefficients are determined such that This ensures the approximate unbiasedness of in estimating . The coefficients can be obtained by equating the coefficients of the corresponding terms of the polynomials on the two sides of (33). To do that, we used the computer algebraic system Maple 16.

Although this computer algebraic system does not give explicit symbolic formulas for the raw moments of a non-central distribution without using special functions, it can efficiently give simple explicit forms of the raw moments of a non-central chi-square distribution up to any required order . The former can then be obtained from the latter by using the following straightforward relationship between the corresponding raw moments of the two distributions (Bain, 1969): where is the th raw moment.

It has to be noted here that the th raw moment of a non-central distribution is finite only for . This puts a constraint on the valid number of polynomial terms to be used in the proposed approximation. If is the size of the control sample, then must be strictly less than .

Now, to apply the approach outlined in Eqs. (31), (32), (33), it remains to find a suitable polynomial approximation to be used in (31). We use two different approximations, the first is based on Bernstein polynomials while the second is a quadrature polynomial approximation.

Bernstein polynomials approximation

From Weierstrass’s Theorem, any continuous real valued function defined on a closed interval can be approximated by a polynomial function. See, for example, Lorentz (1986). In 1912, Bernstein gave a simple probabilistic constructive proof for Weierstrass’s Theorem by introducing the Bernstein polynomials as a series of polynomials that converge uniformly to any continuous bounded function on the closed interval as . See, for example, Chapter (7) in Phillips (2003). The th Bernstein polynomial for is defined as: The polynomial function uniformly approximates on in the sense that (e.g. Theorem 1, Section VII.2 in Feller, 1965).

We use Bernstein polynomials to obtain a polynomial approximation for the chi-square cdf to be used in (31). Although the domain of the chi-square cdf is , we use an affine transformation , for any two arbitrary values and , so as to work on the [0, 1] interval. The two end-points and are initially chosen such that the probability of getting a sample value of the non-centrality parameter outside the interval is fairly negligible. Therefore, we initially take and where, as before, is the value of for which is the -quantile of . As will be shown at the end of Section  3.4, the accuracy of the polynomial approximation is influenced by the choice of and . For extremely large values of , say above the 0.9999-quantile of the chi-square distribution, polynomial functions of small degree are not guaranteed to give a good approximation. Also, accuracy is greatly enhanced if is chosen to be just below the sample median value of . Hence, as a rule of thumb, if is greater than the mode of the chi-square distribution, our final choice of is .

We then approximate in (31) by its th Bernstein polynomial in of the form Clearly, the above expression of is a polynomial of degree in , and we denote its coefficients by . The explicit form of these coefficients was obtained using the computer algebraic system Maple 16. The coefficients of on the left hand side of (33) are equated to their corresponding coefficients in the Bernstein polynomial approximation (36) so as to obtain the values of and, hence, in Eq. (32). In this paper, we obtain the estimate for , 7 and 10, and denote it by , and , respectively.

Quadrature polynomial approximation

We adopt the quadrature formula of Sahai et al. (2004) to obtain another polynomial approximation for . This quadrature formula gives polynomial approximations to the integration of real valued continuous functions defined on the closed interval [0, 1]. Specifically, a function on [0, 1] is approximated by a polynomial in of degree as where partition the interval [0, 1] into equal segments. The polynomial function can then be easily integrated over any sub-interval in [0, 1]. It has been shown empirically that the approximation has good accuracy even when is small. For example, it was used by Richard et al. (2010) to obtain an efficient polynomial approximation of degree 9 to the normal distribution function and its inverse function.

Here, we apply the quadrature formula to approximate the density function of a chi-square distribution on degrees of freedom with a polynomial of degree , which then yields the approximation in (31) after a straightforward symbolic integration of the polynomial.

To work on the [0, 1] interval, we adopt the approach discussed in Section  3.4.1 for choosing the two end-points and , with the same affine transformation . The probability can now be approximated as where is the density function of a chi-square random variable with degrees of freedom.

The coefficients of in the expression of in (38) above are again denoted by , with their explicit forms obtained using Maple 16. The coefficients of in the left hand side of (33) are equated to their corresponding coefficients in the quadrature polynomial approximation in (38) so as to obtain . This gives another form of (Eq. (32)). Here we determine it for , 7 and 10, and denote the resulting estimators by , and , respectively. Simulation results show that these estimators are usually marginally better than those based on Bernstein polynomials.

In general, polynomial functions give approximations that are accurate only on specific intervals of the domain of the underlying approximated function. The accuracy of the polynomial approximations that we use is highly related to the values selected as the two interval end-points, and . Fig. 3, Fig. 4, Fig. 5 show Bernstein and quadrature polynomial approximations for three choices of . The cdf of a chi-square distribution on degrees of freedom is plotted together with its Bernstein approximation and quadrature polynomial approximation, each of degree .

Fig. 3

.

Fig. 4

.

Fig. 5

.

Fig. 3 shows that both polynomial approximations give good accuracy when is the rather short interval [0, 18] of the cdf domain. The value is near the boundary of plausible values for a variate, as 18 is the 0.999 quantile of a distribution. For extremely large values of , neither polynomial approximation of degree 7 is expected to attain good accuracy. This can be seen in Fig. 4, where . For the same extreme value of , if is above the mode of this chi-square distribution (i.e.  ), remarkably better accuracy is obtained, especially for the quadrature approximation. This is shown in Fig. 5, where and . This argument motivates our choice of and that was discussed in Section  3.4.1, as the behavior of both polynomial approximations tends to be the same for all values of . As seen in Fig. 3 and Fig. 5, the Bernstein polynomial approximation can be more accurate than the quadrature polynomial approximation for values near , but the quadrature approximation overall attains better accuracy over the whole interval .

Simulation results

In practice, we expect , the number of test scores on each individual, to be small, while , the number of people in the control sample, may be large. Therefore, we conducted a simulation study that examined combinations of , 4 and 8, with , 20 and 80. Results of some other combinations are available on request from the authors. Based on samples for each combination, we tested the performance of each of the proposed estimators in terms of their average bias, , and the root of mean square error, , denoted by SE (standard error) in tables.

Mahalanobis distance is invariant under affine transformations of location and scale parameters. Since all the methods proposed in this paper depend only on the sample Mahalanobis index, , we chose, without loss of generality, to set the population mean () equal to and the population variance () equal to the identity matrix . The true values of that we examined were 1%, 2.5%, 5%, 10%, 20% and 40%. These true probabilities were attained from Eqs. (2), (5) by choosing each corresponding case’s profile of scores, , as a vector of equal elements. We also examined some cases where the scores in the profile, , are not necessarily equal to each other. But these cases gave almost identical results to those of the profiles with equal scores. We therefore give simulation results only for profiles with equal scores.

Table 1 shows the simulation results for and 4, and . When , , and are the best three estimators in terms of the bias and root mean square error. is slightly better than and up to , but for and 40%, the estimators and are remarkably better than in terms of their bias, with being the best. At , the table shows that , and are again the best three estimators, with showing less bias than for true values of of 10% or more. This suggests that, at small values of , the best two competitor estimates are and , where the former is doing better at smaller values of the true probability .

Table 1

Bias and root mean square error (SE) of the proposed estimates of at .

ν_1		1%		2.5%		5%		10%		20%		40%
		Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE
2	P^_F	4.6	7.1	5.9	9.4	6.8	11.5	7.2	14.0	6.2	16.5	2.3	18.7
	P^_χ^2	0.3	2.5	0.1	4.1	−0.5	6.1	−1.8	9.3	−4.5	14.2	−9.0	21.2
	P^_D	1.8	4.6	2.5	6.9	3.0	9.3	3.3	12.5	2.7	16.5	0.3	20.6
	P^_MD	1.8	4.6	2.5	6.9	3.0	9.3	3.3	12.4	2.4	16.0	−1.0	19.1
	P^_BY	4.5	6.9	5.5	8.9	6.4	11.1	7.2	13.9	5.4	16.1	2.1	19.1
	P^_M	3.4	6.8	4.9	9.7	6.4	12.7	7.8	16.3	8.6	20.1	7.5	22.8
	P^_R	6.8	10.6	9.4	14.2	11.8	17.5	14.0	20.9	15.1	23.5	12.5	22.8
	P^_T	0.8	3.9	1.2	6.2	1.6	9.1	2.0	13.0	2.1	18.1	1.7	22.9
	P^_B4	2.0	5.0	2.8	7.4	3.6	10.1	4.2	13.7	4.3	18.1	3.2	22.1
	P^_Q4	0.9	3.8	1.3	6.2	1.6	8.9	1.9	12.8	1.9	17.8	1.3	22.8
4	P^_F	7.0	10.5	8.6	13.2	9.7	15.4	10.1	17.6	8.4	19.4	2.5	20.6
	P^_χ^2	−0.2	2.2	−0.9	3.7	−2.1	5.7	−4.7	9.3	−9.6	15.7	−18.1	26.2
	P^_D	2.8	7.1	3.9	9.9	4.7	12.8	5.1	16.2	4.1	20.2	0.3	23.9
	P^_MD	2.8	7.0	3.9	9.8	4.7	12.5	4.9	15.7	3.4	18.7	−2.1	21.1
	P^_BY	7.2	10.7	8.4	13.2	9.7	15.3	9.4	17.4	8.2	19.4	2.9	20.5
	P^_M	5.3	10.4	7.4	14.1	9.2	17.5	10.8	21.2	11.2	24.7	8.4	26.1
	P^_R	10.6	16.4	14.0	20.9	17.0	24.7	19.7	28.3	20.6	30.5	16.9	28.2
	P^_T	1.6	6.4	2.4	9.7	3.3	13.4	4.1	18.0	4.5	23.7	3.5	28.4
	P^_B4	3.1	7.6	4.3	10.7	5.4	13.8	6.1	17.7	5.9	22.1	3.1	25.5
	P^_Q4	1.4	5.8	2.0	8.8	2.6	12.2	2.9	16.6	2.7	22.5	1.7	28.7

An important point from Table 1 is the cautionary message that each of the methods shows noticeable bias for some values of at some combinations of and . Four methods perform particularly poorly: , , , .

Table 2 shows the simulation results at , 4 and 8, and . For all listed values of , the three estimates , and are doing better than the others. For small values of , the bias of both and is generally less than that of , while the root mean square error of is always less than those of and . Comparing to , it can be seen in Table 2 that is better than in terms of bias, although the root mean square error of is slightly greater than that of for all values of the true probability except . This suggests that is the best estimate at as it has rather small values of both bias and root mean square error for all values of the true probability .

Table 2

Bias and root mean square error (SE) of the proposed estimates of at .

ν_1		1%		2.5%		5%		10%		20%		40%
		Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE
2	P^_F	2.2	3.8	3.1	5.6	3.7	7.4	4.1	9.6	3.5	12.2	1.2	14.6
	P^_χ^2	0.2	1.8	0.1	3.1	−0.3	4.8	−1.1	7.4	−2.6	11.2	−5.1	15.9
	P^_D	0.9	2.6	1.3	4.3	1.6	6.2	1.8	8.8	1.3	12.1	0.0	15.5
	P^_MD	0.9	2.6	1.4	4.3	1.7	6.2	1.8	8.8	1.4	12.1	−0.3	15.0
	P^_BY	2.2	3.8	3.2	5.8	3.8	7.3	4.1	9.4	3.4	12.0	1.5	14.7
	P^_M	1.5	3.3	2.3	5.3	3.0	7.4	3.8	10.2	4.2	13.5	3.6	16.3
	P^_R	2.4	4.4	3.7	6.6	4.9	8.9	6.2	11.8	6.8	14.5	5.6	16.0
	P^_T	0.3	2.1	0.4	3.8	0.5	5.8	0.6	8.8	0.5	12.8	0.4	16.5
	P^_B4	0.9	2.6	1.3	4.4	1.7	6.4	2.0	9.2	1.9	12.7	1.4	16.2
	P^_Q4	0.2	2.0	0.2	3.7	0.3	5.8	0.4	8.8	0.3	12.7	0.2	16.6
	P^_B7	0.6	2.4	0.9	4.1	1.2	6.1	1.4	9.0	1.3	12.7	0.9	16.3
	P^_Q7	−0.1	1.9	−0.1	3.6	−0.1	5.7	0.0	8.8	−0.1	12.9	0.0	16.7
4	P^_F	3.5	5.7	4.7	8.0	5.6	10.1	6.0	12.4	5.1	14.8	1.3	16.8
	P^_χ^2	−0.1	1.8	−0.5	3.2	−1.3	5.1	−3.1	8.1	−6.3	13.1	−11.4	20.3
	P^_D	1.4	3.9	2.1	6.0	2.6	8.3	2.8	11.4	2.3	15.0	−0.1	18.4
	P^_MD	1.4	3.9	2.1	6.1	2.7	8.4	2.9	11.3	2.1	14.7	−1.0	17.3
	P^_BY	3.2	5.4	4.8	7.9	5.7	9.8	5.4	11.6	4.5	15.2	2.2	16.8
	P^_M	2.2	4.9	3.4	7.5	4.5	10.1	5.4	13.4	5.7	16.9	4.1	19.2
	P^_R	3.6	6.5	5.4	9.5	7.0	12.4	8.6	15.7	9.4	18.8	7.5	19.5
	P^_T	0.5	3.1	0.8	5.3	1.0	7.9	1.3	11.6	1.3	16.3	0.8	20.5
	P^_B4	1.4	3.8	2.0	6.1	2.6	8.5	3.0	11.8	2.8	15.8	1.2	19.3
	P^_Q4	0.2	2.8	0.4	5.0	0.5	7.6	0.6	11.3	0.5	16.1	0.0	20.7
	P^_B7	0.9	3.4	1.4	5.6	1.8	8.1	2.1	11.6	2.0	15.9	0.8	19.8
	P^_Q7	−0.1	2.6	−0.1	4.9	0.0	7.6	0.0	11.5	0.1	16.5	0.0	21.0
8	P^_F	5.6	8.9	7.3	11.6	8.5	14.0	9.0	16.4	7.5	18.4	2.1	19.7
	P^_χ^2	−0.6	1.4	−1.5	2.8	−3.1	4.9	−6.3	8.7	−12.3	15.7	−22.2	27.4
	P^_D	2.4	6.1	3.4	8.9	4.2	11.8	4.6	15.2	3.8	19.1	0.4	22.6
	P^_MD	2.4	6.1	3.4	8.9	4.2	11.7	4.5	14.8	3.2	17.8	−1.5	20.2
	P^_BY	5.7	9.3	7.4	11.7	8.2	13.8	8.4	16.0	7.1	18.4	2.8	19.2
	P^_M	3.6	7.9	5.3	11.2	6.8	14.5	8.1	18.1	8.2	21.7	5.5	23.7
	P^_R	5.7	10.5	8.2	14.4	10.4	17.9	12.5	21.7	13.3	24.7	10.6	24.8
	P^_T	1.0	5.1	1.6	8.2	2.3	11.6	2.8	16.1	3.0	21.6	2.0	26.3
	P^_B4	2.2	6.0	3.2	8.9	4.1	11.9	4.6	15.6	4.2	19.9	1.3	23.4
	P^_Q4	0.5	4.4	0.8	7.3	1.1	10.6	1.2	15.0	0.8	20.6	−0.7	25.6
	P^_B7	1.6	5.4	2.3	8.3	3.0	11.4	3.4	15.4	3.0	20.3	0.9	24.6
	P^_Q7	0.1	4.2	0.2	7.2	0.4	10.8	0.5	15.6	0.6	22.0	0.4	27.7

Simulation results for (, 4 and 8) are presented in Table 3. This value of is large enough for the estimators and to be computed. All estimators are doing well at this very large sample size and all have similar bias and root mean square error. However, the estimators based on polynomial approximations are still slightly better than the others. Specifically, the quadrature based estimators , and have very low bias as shown in Table 3, with the bias of always less than or equal to those of and .

Table 3

Bias and root mean square error (SE) of the proposed estimates of at .

ν_1		1%		2.5%		5%		10%		20%		40%
		Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE	Bias	SE
2	P^_F	0.5	1.2	0.8	2.0	1.0	3.1	1.2	4.5	1.0	6.4	0.3	8.1
	P^_χ^2	0.1	0.8	0.0	1.6	0.0	2.6	−0.3	4.2	−0.8	6.2	−1.4	8.3
	P^_D	0.2	1.0	0.3	1.8	0.5	2.8	0.5	4.4	0.4	6.4	0.0	8.2
	P^_MD	0.2	1.0	0.4	1.8	0.4	2.9	0.5	4.4	0.4	6.4	0.0	8.2
	P^_BY	0.6	1.2	0.8	2.0	1.1	3.2	1.2	4.3	0.9	6.4	0.2	8.1
	P^_M	0.3	1.1	0.5	1.9	0.8	3.0	1.0	4.6	1.1	6.6	0.9	8.3
	P^_R	0.5	1.2	0.8	2.1	1.1	3.2	1.4	4.8	1.6	6.7	1.3	8.2
	P^_T	0.0	0.9	0.0	1.7	0.0	2.8	0.1	4.4	0.0	6.5	0.0	8.4
	P^_B4	0.2	0.9	0.3	1.8	0.4	2.8	0.4	4.4	0.4	6.5	0.3	8.3
	P^_Q4	0.0	0.8	0.0	1.7	0.0	2.8	0.0	4.4	0.0	6.5	0.0	8.4
	P^_B7	0.1	0.9	0.2	1.7	0.2	2.8	0.3	4.4	0.3	6.5	0.2	8.3
	P^_Q7	0.0	0.8	0.0	1.7	0.0	2.8	0.0	4.4	0.0	6.5	0.0	8.4
	P^_B10	0.1	0.9	0.1	1.7	0.2	2.8	0.2	4.4	0.2	6.5	0.1	8.3
	P^_Q10	0.0	0.8	0.0	1.7	0.0	2.8	0.0	4.4	0.0	6.5	0.0	8.4
4	P^_F	0.8	1.7	1.3	2.8	1.6	4.1	1.8	5.8	1.5	8.0	0.4	9.8
	P^_χ^2	0.0	1.0	−0.1	1.9	−0.4	3.1	−1.0	5.0	−2.0	7.7	−3.5	10.5
	P^_D	0.4	1.3	0.6	2.4	0.7	3.7	0.8	5.5	0.6	8.0	−0.1	10.1
	P^_MD	0.4	1.3	0.6	2.4	0.7	3.7	0.8	5.6	0.6	8.0	−0.1	9.9
	P^_BY	0.8	1.7	1.2	2.8	1.6	4.1	1.6	5.8	1.3	7.8	0.4	9.9
	P^_M	0.5	1.4	0.8	2.6	1.1	3.9	1.4	5.8	1.5	8.3	1.0	10.2
	P^_R	0.7	1.6	1.1	2.8	1.5	4.2	2.0	6.1	2.2	8.5	1.7	10.2
	P^_T	0.0	1.1	0.1	2.1	0.1	3.5	0.1	5.5	0.1	8.2	0.0	10.4
	P^_B4	0.3	1.2	0.4	2.3	0.6	3.6	0.6	5.6	0.6	8.1	0.3	10.3
	P^_Q4	0.0	1.1	0.0	2.1	0.0	3.5	0.0	5.5	0.0	8.2	0.0	10.4
	P^_B7	0.2	1.2	0.3	2.2	0.4	3.6	0.4	5.5	0.4	8.1	0.1	10.3
	P^_Q7	0.0	1.1	0.0	2.1	0.0	3.5	0.0	5.5	0.0	8.2	0.0	10.4
	P^_B10	0.1	1.1	0.2	2.2	0.3	3.5	0.3	5.5	0.3	8.2	0.1	10.4
	P^_Q10	0.0	1.1	0.0	2.1	0.0	3.5	0.0	5.5	0.0	8.2	0.0	10.4
8	P^_F	1.4	2.6	2.0	4.1	2.6	5.7	2.9	7.8	2.6	10.3	0.6	12.3
	P^_χ^2	−0.2	1.1	−0.6	2.1	−1.3	3.6	−2.6	6.0	−4.9	9.9	−8.1	14.6
	P^_D	0.6	1.9	0.9	3.3	1.2	5.0	1.4	7.3	1.1	10.3	0.0	13.0
	P^_MD	0.6	1.9	0.9	3.3	1.2	5.0	1.4	7.3	1.1	10.3	−0.3	12.5
	P^_BY	1.4	2.6	1.9	3.9	2.3	5.6	2.8	7.9	2.3	10.1	0.9	12.2
	P^_M	0.8	2.1	1.3	3.6	1.7	5.4	2.1	7.8	2.2	10.7	1.4	13.1
	P^_R	1.0	2.3	1.7	4.0	2.3	5.8	3.0	8.3	3.3	11.1	2.5	13.2
	P^_T	0.1	1.5	0.2	2.9	0.2	4.6	0.3	7.2	0.3	10.7	0.1	13.7
	P^_B4	0.5	1.8	0.7	3.2	0.9	4.9	1.1	7.3	0.9	10.5	0.3	13.3
	P^_Q4	0.0	1.5	0.0	2.8	0.0	4.6	0.0	7.2	0.0	10.7	0.0	13.7
	P^_B7	0.3	1.6	0.5	3.0	0.6	4.7	0.7	7.2	0.6	10.5	0.1	13.5
	P^_Q7	0.0	1.4	0.0	2.8	0.0	4.6	0.0	7.2	0.0	10.7	0.0	13.8
	P^_B10	0.2	1.6	0.4	3.0	0.5	4.7	0.5	7.2	0.4	10.6	0.1	13.6
	P^_Q10	0.0	1.4	0.0	2.8	0.0	4.6	0.0	7.2	0.0	10.7	0.0	13.8

Feasibility of estimates and absolute error

Ranges of estimates

The simulations in Section  4 show that the estimators based on polynomial approximations performed well, in that they had the minimum bias among the reported estimators. However, this does not mean that the estimates they produce are always sensible. Specifically, when is very small they can give estimates of the proportion that are greater than 100%, and when is very big they can give estimates that are less than zero. This problem comes to light by studying the behavior of the proposed estimators at different values in the domain of .

For example, at and , Fig. 6(a) shows that the twelve estimators all have similar patterns for . But a closer look at the part of the domain where (Fig. 6(b)) reveals that , , and are not monotonically decreasing with and they exceed 100% at some values of . With the same degrees of freedom, another problem appears in Fig. 6(c), where both and are below zero for some values of .

Fig. 6

Estimates of given by different methods for various combinations of and .

Similar problems appear at large sample sizes as well. For example, at and , is slightly below zero for some values of and as shown in Fig. 6(d), , and all exceed 100% for some values of . The problem does not arise with other estimators—estimates are always in the range 0%–100% for , , , , , , and . However, as Fig. 6(b) shows, the estimate of approaches 100% as approaches 0. This is clearly unrealistic as the case’s value will not equal the population mean , even if equals the sample mean . As in Section  3.1.2, a pragmatic approach is to treat the case’s profile as that of a randomly chosen control when the case’s profile seems nearer to than would be expected of a control’s profile.

Performances as measured by absolute error

In Section  4, mean square error and average error (bias) were used to evaluate the performance of the various estimators considered in this paper. Alternative evaluation criteria include average absolute error () and median error (ME  =  median()). Here, we briefly examine the performance of our estimators under these criteria. In theory, the median estimator should give a median error of 0 and have a lower AAE than other estimators whose median error is small. Its closely related estimator, , should also perform well.

Results for and are presented in Table 4. It can be seen that the median error in estimating is 0.0 for both and for each of the tabulated values of . In contrast, the median error of every other estimator is never 0.0 except for when ; otherwise the median error of the other estimators is typically quite marked.

Table 4

Median error (ME) and average absolute error (AAE) of the median estimates of at and .

	1%		2.5%		5%		40%
	ME	AAE	ME	AAE	ME	AAE	ME	AAE
P^_F	1.6	3.0	2.3	4.6	2.9	6.2	1.1	13.3
P^_χ^2	−0.8	0.9	−1.9	2.1	−3.6	3.7	−16.0	17.8
P^_D	0.0	1.8	0.0	3.3	0.0	5.1	0.0	14.5
P^_MD	0.0	1.8	0.0	3.3	0.0	5.1	0.0	13.7
P^_BY	1.4	2.8	2.2	4.8	2.9	6.5	1.1	13.2
P^_M	0.3	2.2	0.6	3.9	1.0	5.8	3.0	14.9
P^_R	0.8	2.7	1.4	4.6	2.3	6.7	5.6	15.1
P^_T	−0.7	1.4	−1.3	2.9	−2.0	4.9	0.4	16.1
P^_B4	−0.1	1.7	−0.3	3.2	−0.4	5.0	0.4	15.1
P^_Q4	−0.9	1.3	−1.7	2.9	−2.5	4.8	−0.3	16.2
P^_B7	−0.4	1.5	−0.8	3.0	−1.1	4.9	0.1	15.5
P^_Q7	−1.0	1.3	−1.9	2.8	−2.8	4.9	−0.1	16.4
P^_B10	−0.5	1.5	−1.1	2.9	−1.5	4.9	0.0	15.8
P^_Q10	−0.9	1.3	−1.8	2.8	−2.8	4.8	−0.1	16.4

A fuller examination of the median errors given by and is provided in Table 5, where results are given for these estimators for all the combinations of and that were considered in Table 1, Table 2, Table 3. The two estimators give identical median error for every combination and that error is very small in every case. Hence, if we want an estimator that has very small median error, then both and can fill that role. The average absolute error is marginally better with the estimator, but the differences are very slight. However, is the preferable estimator because it will not give unrealistic estimates of , while will sometimes estimate as 100% when that is not a credible estimate. Consequently, if a point estimator of is required, one reasonable choice is to give as the estimator and say that it gives small median error without making any claim about its bias (average error).

Table 5

Median error (ME) and average absolute error (AAE) of the median estimates and of .

ν_2	ν_1		1%		2.5%		5%		40%
			ME	AAE	ME	AAE	ME	AAE	ME	AAE
10	2	P^_D	0.0	2.5	0.0	4.2	0.0	6.3	−0.1	17.1
	P^_MD	0.0	2.5	0.0	4.2	0.0	6.3	−0.1	16.1
	4	P^_D	0.0	3.6	0.0	5.8	0.0	8.3	0.0	20.2
	P^_MD	0.0	3.6	0.0	5.8	0.0	8.2	0.0	17.8
20	2	P^_D	0.0	1.6	0.0	2.8	0.0	4.4	0.0	12.6
	P^_MD	0.0	1.6	0.0	2.8	0.0	4.4	0.0	12.3
	4	P^_D	0.0	2.1	0.0	3.7	0.0	5.7	0.0	15.3
	P^_MD	0.0	2.1	0.0	3.7	0.0	5.7	0.0	14.3
	8	P^_D	0.0	3.1	0.0	5.2	−0.1	7.6	−0.1	19.1
	P^_MD	0.0	3.1	0.0	5.2	−0.1	7.5	−0.1	17.1
80	2	P^_D	0.0	0.7	0.0	1.3	0.0	2.2	0.1	6.6
	P^_MD	0.0	0.7	0.0	1.3	0.0	2.2	0.1	6.6
	4	P^_D	0.0	0.9	0.0	1.7	0.0	2.8	0.0	8.2
	P^_MD	0.0	0.9	0.0	1.7	0.0	2.8	0.0	8.1
	8	P^_D	0.0	1.2	0.0	2.2	0.0	3.6	0.0	10.6
	P^_MD	0.0	1.2	0.0	2.2	0.0	3.6	0.0	10.2

Concluding comments

The task that motivated this paper seemed straightforward: find a good point estimator of the abnormality of a Mahalanobis index. The answer is less straightforward, as the best choice of estimator will depend on the purpose for which the estimator is required. The following summarizes our findings.

The most common criteria used to choose an estimator are bias and mean square error; the minimum variance unbiased estimator is often the preferred estimator if such an estimator can be found. Under these criteria the best estimators are those based on a quadrature polynomial approximation, , and , provided occasional negative estimates are not a problem. (The negative estimates would presumably be set to 0.) Only can be used for ; is best for ; and are marginally the best ( is almost as good) for .

If mean square error is to be minimized and bias is unimportant, then is the best estimator, but it displays substantial bias even when is large.

Sometimes, an estimate of is to be used as an input into further analysis. Commonly though, an estimate of is to be communicated to others (perhaps in a journal paper or a technical report) and then a good descriptive statistic is required. In that context, the best estimator would seem to be the modified median estimator, . It should be referred to as the median estimator as that is accurate: it is designed to give low median bias rather than low average bias and, indeed, its median bias is very low. It is preferable to the median estimate () because it always gives sensible estimates while sometimes gives estimates that are unrealistically small when judged by common sense.

Based on our simulation results, we recommend that should generally be used as the point estimator of . However, if unbiasedness of the required estimate is crucially important we recommend that should be used for and should be used for . Out-of-range values of these two estimators need to be artificially constrained so as not to lie outside the interval [0, 1].

This work is part of an on-going project that develops statistical methods for analyzing single patient data, and these recommendations are implemented in software for making inferences from Mahalanobis distance about the abnormality of an individual’s test score profile. Previous methods that we have developed are well-used by neuropsychologists (see, for example, papers that cite Crawford and Garthwaite, 2002, Crawford and Garthwaite, 2005, Crawford and Garthwaite, 2007) so it is likely that the recommendations will influence practice. The work in this paper makes these recommendations well-informed.