ExploringDataBlog

Measuring associations between non-numeric variables

2012-02-04T16:06:00.000-08:00

It is often useful to know how strongly or weakly two variables are associated: do they vary together or are they essentially unrelated? In the case of numerical variables, the best-known measure of association is the product-moment correlation coefficient introduced by Karl Pearson at the end of the nineteenth century. For variables that are ordered but not necessarily numeric (e.g., Likert scale responses with levels like “strongly agree,” “agree,” “neither agree nor disagree,” “disagree” and “strongly disagree”), association can be measured in terms of the Spearman rank correlation coefficient. Both of these measures are discussed in detail in Chapter 10 of Exploring Data in Engineering, the Sciences, and Medicine. For unordered categorical variables (e.g., country, state, county, tumor type, literary genre, etc.), neither of these measures are applicable, but applicable alternatives do exist. One of these is Goodman and Kruskal’s tau measure, discussed very briefly in Exploring Data (Chapter 10, page 492). The point of this post is to give a more detailed discussion of this association measure, illustrating some of its advantages, disadvantages, and peculiarities.

A more complete discussion of Goodman and Kruskal’s tau measure is given in Agresti’s book Categorical Data Analysis, on pages 68 and 69. It belongs to a family of categorical association measures of the general form:

a(x,y) = [V(y) – E{V(y|x)}]/V(y)

where V(y) is a measure of the overall (i.e., marginal) variability of y and E{V(y|x)} is the expected value of the conditional variability V(y|x) of y given a fixed value of x, where the expectation is taken over all possible values of x. These variability measures can be defined in different ways, leading to different association measures, including Goodman and Kruskal’s tau as a special case. Agresti’s book gives detailed expressions for several of these variability measures, including the one on which Goodman and Kruskal’s tau is based, and an alternative expression for the overall association measure a(x,y) is given in Eq. (10.178) on page 492 of Exploring Data. This association measure does not appear to be available in any current R package, but it is easily implemented as the following function:

GKtau <- function(x,y){
#
# First, compute the IxJ contingency table between x and y
#
Nij = table(x,y,useNA="ifany")
#
# Next, convert this table into a joint probability estimate
#
PIij = Nij/sum(Nij)
#
# Compute the marginal probability estimates
#
PIiPlus = apply(PIij,MARGIN=1,sum)
PIPlusj = apply(PIij,MARGIN=2,sum)
#
# Compute the marginal variation of y
#
Vy = 1 - sum(PIPlusj^2)
#
# Compute the expected conditional variation of y given x
#
InnerSum = apply(PIij^2,MARGIN=1,sum)
VyBarx = 1 - sum(InnerSum/PIiPlus)
#
# Compute and return Goodman and Kruskal's tau measure
#
tau = (Vy - VyBarx)/Vy
tau
}

An important feature of this procedure is that it allows missing values in either of the variables x or y, treating “missing” as an additional level. In practice, this is sometimes very important since missing values in one variable may be strongly associated with either missing values in another variable or specific non-missing levels of that variable.

An important characteristic of Goodman and Kruskal’s tau measure is its asymmetry: because the variables x and y enter this expression differently, the value of a(y,x) is not the same as the value of a(x,y), in general. This stands in marked contrast to either the product-moment correlation coefficient or the Spearman rank correlation coefficient, which are both symmetric, giving the same association between x and y as that between y and x. The fundamental reason for the asymmetry of the general class of measures defined above is that they quantify the extent to which the variable x is useful in predicting y, which may be very different than the extent to which the variable y is useful in predicting x. Specifically, if x and y are statistically independent, then E{V(y|x)} = V(y) – i.e., knowing x does not help at all in predicting y – and this implies that a(x,y) = 0. At the other extreme, if y is perfectly predictable from x, then E{V(y|x)} = 0, which implies that a(x,y) = 1. As the examples presented next demonstrate, it is possible that y is extremely predictable from x, but x is only slightly predictable from y.

Specifically, consider the sequence of 400 random numbers, uniformly distributed between 0 and 1 generated by the following R code:

set.seed(123)

u = runif(400)

(Here, I have used the “set.seed” command to initialize the random number generator so repeated runs of this example will give exactly the same results.) The second sequence is obtained by quantizing the first, rounding the values of u to a single digit:

x = round(u,digits=1)

The plot below shows the effects of this coarse quantization: values of u vary continuously from 0 to 1, but values of x are restricted to 0.0, 0.1, 0.2, … , 1.0. Although this example is simulation-based, it is important to note that this type of grouping of variables is often encountered in practice (e.g., the use of age groups instead of ages in demographic characterizations, blood pressure characterizations like “normal,” “borderline hypertensive,” etc. in clinical data analysis, or the recording of industrial process temperatures to the nearest 0.1 degree, in part due to measurement accuracy considerations and in part due to memory limitations of early data collection systems).

In this particular case, because the variables x and u are both numeric, we could compute either the product-moment correlation coefficient or the Spearman rank correlation, obtaining the very large value of approximately 0.995 for either one, showing that these variables are strongly associated. We can also apply Goodman and Kruskal’s tau measure here, and the result is much more informative. Specifically, the value of a(u,x) is 1 in this case, correctly reflecting the fact that the grouped variable x is exactly computable from the original variable u. In contrast, the value of a(x,u) is approximately 0.025, suggesting – again correctly – that the original variable u cannot be well predicted from the grouped variable x.

To illustrate a case where the product-moment and rank correlation measures are not applicable at all, consider the following alphabetic re-coding of the variable x into an unordered categorical variable c:

letters = c(“A”, “B”, “C”, “D”, “E”, “F”, “G”, “H”, “I”, “J”, “K”)

c = letters[10*x+1]

In this case, both of the Goodman and Kruskal tau measures, a(x,c) and a(c,x), are equal to 1, reflecting the fact that these two variables are effectively identical, related via the non-numeric transformation given above.

Being able to detect relationships like these can be extremely useful in exploratory data analysis where such relationships may be unexpected, particularly in the early stages of characterizing a dataset whose metadata – i.e., detailed descriptions of the variables included in the dataset – is absent, incomplete, ambiguous, or suspect. As a real data illustration, consider the rent data frame from the R package gamlss.data, which has 1,969 rows, each corresponding to a rental property in Munich, and 9 columns, each giving a characteristic of that unit (e.g., the rent, floor space, year of construction, etc.). Three of these variables are Sp, a binary variable indicating whether the location is considered above average (1) or not (0), Sm, another binary variable indicating whether the location is considered below average (1) or not (0), and loc, a three-level variable combining the information in these other two, taking the values 1 (below average), 2 (average), or 3 (above average). The Goodman and Kruskal tau values between all possible pairs of these three variables are:

a(Sm,Sp) = a(Sp,Sm) = 0.037

a(Sm,loc) = 0.245 vs. a(loc,Sm) = 1

a(Sp,loc) = 0.701 vs. a(loc,Sp) = 1

The first of these results – the symmetry of Goodman and Kruskal’s tau for the variables Sm and Sp – is a consequence of the fact that this measure is symmetric for any pair of binary variables. In fact, the odds ratio that I have discussed in previous posts represents a much better way of characterizing the relationship between binary variables (here, the odds ratio between Sm and Sp is zero, reflecting the fact that a location cannot be both “above average” and “below average” at the same time). The real utility of the tau measure here is that the second and third lines above show that the variables Sm and Sp are both re-groupings of the finer-grained variable loc.

Finally, a more interesting exploratory application to this dataset is the following one. Computing Goodman and Kruskal’s tau measure between the location variable loc and all of the other variables in the dataset – beyond the cases of Sm and Sp just considered – generally yields small values for the associations in either direction. As a specific example, the association a(loc,Fl) is 0.001, suggesting that location is not a good predictor of the unit’s floor space in meters, and although the reverse association a(Fl,loc) is larger (0.057), it is not large enough to suggest that the unit’s floor space is a particularly good predictor of its location quality. The same is true of most of the other variables in the dataset: they are neither well predicted by nor good predictors of location quality. The one glaring exception is the rent variable R: although the association a(loc,R) is only 0.001, the reverse association a(R,loc) is 0.907, a very large value suggesting that location quality is quite well predicted by the rent. The beanplot above shows what is happening here: because the variation in rents for all three location qualities is substantial, knowledge of the loc value is not sufficient to accurately predict the rent R, but these rent values do generally increase in going from below-average locations (loc = 1) to average locations (loc = 2) to above-average locations (loc = 3). For comparison, the beanplots below show why the association with floor space is so much weaker: both the mean floor space in each location quality group and the overall range of these values are quite comparable, implying that neither location quality can be well predicted from floor space nor vice versa.

The asymmetry of Goodman and Kruskal’s tau measure is disconcerting at first because it has no counterpart in better-known measures like the product-moment correlation coefficient between numerical variables, Spearman’s rank correlation coefficient between ordinal variables, or the odds ratio between binary variables. One of the points of this post has been to demonstrate how this unusual asymmetry can be useful in practice, distinguishing between the ability of one variable x to predict another variable y, and the reverse case.

Moving window filters and the pracma package

2012-01-14T11:06:00.000-08:00

In my last post, I discussed the Hampel filter, a useful moving window nonlinear data cleaning filter that is available in the R package pracma. In this post, I briefly discuss this moving window filter in a little more detail, focusing on two important practical points: the choice of the filter’s local outlier detection threshold, and the question of how to initialize moving window filters. This second point is particularly important here because the pracma package initializes the Hampel filter in a particularly appropriate way, but doesn’t do such a good job of initializing the Savitzky-Golay filter, a linear smoothing filter that is popular in physics and chemistry. Fortunately, this second difficulty is easy to fix, as I demonstrate here.

Recall from my last post that the Hampel filter is a moving window implementation of the Hampel identifier, discussed in Chapter 7 of Exploring Data in Engineering, the Sciences, and Medicine. In particular, this procedure – implemented as outlierMAD in the pracma package – is a nonlinear data cleaning filter that looks for local outliers in a time-series or other streaming data sequence, replacing them with a more reasonable alternative value when it finds them. Specifically, this filter may be viewed as a more effective alternative to a “local three-sigma edit rule” that would replace any data point lying more than three standard deviations from the mean of its neighbors with that mean value. The difficulty with this simple strategy is that both the mean and especially the standard deviation are badly distorted by the presence of outliers in the data, causing this data cleaning procedure to often fail completely in practice. The Hampel filter instead uses the median of neighboring observations as a reference value, and the MAD scale estimator as an alternative measure of distance: that is, a data point is declared an outlier and replaced if it lies more than some number t of MAD scale estimates from the median of its neighbors; the replacement value used in this procedure is the median.

More specifically, for each observation in the original data sequence, the Hampel filter constructs a moving window that includes the K prior points, the data point of primary interest, and the K subsequent data points. The reference value used for the central data point is the median of these 2K+1 successive observations, and the MAD scale estimate is computed from these same observations to serve as a measure of the “natural local spread” of the data sequence. If the central data point lies more than t MAD scale estimate values from the median, it is replaced with the median; otherwise, it is left unchanged. To illustrate the performance of this filter, the top plot above shows the sequence of 1024 successive physical property measurements from an industrial manufacturing process that I also discussed in my last post. The bottom plot in this pair shows the results of applying the Hampel filter with a window half-width parameter K=5 and a threshold value of t = 3 to this data sequence. Comparing these two plots, it is clear that the Hampel filter has removed the glaring outlier – the value zero – at observation k = 291, yielding a cleaned data sequence that varies over a much narrower (and, at least in this case, much more reasonable) range of possible values. What is less obvious is that this filter has also replaced 18 other data points with their local median reference values.

The above plot shows the original data sequence, but on approximately the same range as the cleaned data sequence so that the glaring outlier at k = 291 no longer dominates the figure. The large solid circles represent the 18 additional points that the Hampel filter has declared to be outliers and replaced with their local median values. This plot was generated using the Hampel filter implemented in the outlierMAD command in the pracma package, which has the following syntax:

outlierMAD(x,k)

where x is the data sequence to be cleaned and k is the half-width that defines the moving data window on which the filter is based. Here, specifying k = 5 results in an 11-point moving data window. Unfortunately, the threshold parameter t is hard-coded as 3 in this pracma procedure, which has the following code:

outlierMAD <- function (x, k){

n <- length(x)

y <- x

ind <- c()

L <- 1.4826

t0 <- 3

for (i in (k + 1):(n - k)) {

x0 <- median(x[(i - k):(i + k)])

S0 <- L * median(abs(x[(i - k):(i + k)] - x0))

if (abs(x[i] - x0) > t0 * S0) {

y[i] <- x0

ind <- c(ind, i)

}

list(y = y, ind = ind)

}

Note that it is a simple matter to create your own version of this filter, specifying the threshold (here, the variable t0) to have a default value of 3, but allowing the user to modify it in the function call. Specifically, the code would be:

HampelFilter <- function (x, k,t0=3){

n <- length(x)

y <- x

ind <- c()

L <- 1.4826

for (i in (k + 1):(n - k)) {

x0 <- median(x[(i - k):(i + k)])

S0 <- L * median(abs(x[(i - k):(i + k)] - x0))

if (abs(x[i] - x0) > t0 * S0) {

y[i] <- x0

ind <- c(ind, i)

}

list(y = y, ind = ind)

}

The advantage of this modification is that it allows you to explore the influence of varying the threshold parameter. Note that increasing t0 makes the filter more forgiving, allowing more extreme local fluctuations to pass through the filter unmodified, while decreasing t0 makes the filter more aggressive, declaring more points to be local outliers and replacing them with the appropriate local median. In fact, this filter remains well-defined even for t0 = 0, where it reduces to the median filter, popular in nonlinear digital signal processing. John Tukey – the developer or co-developer of many useful things, including the fast Fourier transform (FFT) – introduced the median filter at a technical conference in 1974, and it has profoundly influenced subsequent developments in nonlinear digital filtering. It may be viewed as the most aggressive limit of the Hampel filter and, although it is quite effective in removing local outliers, it is often too aggressive in practice, introducing significant distortions into the original data sequence. This point may be seen in the plot below, which shows the results of applying the median filter (i.e., the HampelFilter procedure defined above with t0=0) to the physical property dataset. In particular, the heavy solid line in this plot shows the behavior of the first 250 points of the median filtered sequence, while the lighter dotted line shows the corresponding results for the Hampel filter with t0=3. Note the “clipped” or “blocky” appearance of the median filtered results, compared with the more irregular local variation seen in the Hampel filtered results. In many applications (e.g., fitting time-series models), the less aggressive Hampel filter gives much better overall results.

The other main issue I wanted to discuss in this post is that of initializing moving window filters. The basic structure of these filters – whether they are nonlinear types like the Hampel and median filters discussed above, or linear types like the Savitzky-Golay filter discussed briefly below – is built on a moving data window that includes a central point of interest, prior observations and subsequent observations. For a symmetric window that includes K prior and K subsequent observations, this window is not well defined for the first K or the last K observations in the data sequence. These points must be given special treatment, and a very common approach in the digital signal processing community is to extend the original sequence by appending K additional copies of the first element to the beginning of the sequence and K additional copies of the last element to the end of the sequence. The pracma implementation of the Hampel filter procedure (outlierMAD) takes an alternative approach, one that is particularly appropriate for data cleaning filters. Specifically, procedure outlierMAD simply passes the first and last K observations unmodified from the original data sequence to the filter output. This would also seem to be a reasonable option for smoothing filters like the linear Savitzky-Golay filter discussed next.

As noted, this linear smoothing filter is popular in chemistry and physics, and it is implemented in the pracma package as procedure savgol. For a more detailed discussion of this filter, refer to the treatment in the book Numerical Recipes, which the authors of the pracma package cite for further details (Section 14.8). Here, the key point is that this filter is a linear smoother, implemented as the convolution of the input sequence with an impulse response function (i.e., a smoothing kernel) that is constructed by the savgol procedure. The above two plots show the effects of applying this filter with a total window width of 11 points (i.e., the same half-width K = 5 used with the Hampel and median filters), first to the raw physical property data sequence (upper plot), and then to the sequence after it has been cleaned by the Hampel filter (lower plot). The large downward spike at k = 291 in the upper plot reflects the impact of the glaring outlier in the original data sequence, illustrating the practical importance of removing these artifacts from a data sequence before applying smoothing procedures like the Savitzky-Golay filter. Both the upper and lower plots exhibit similarly large spikes at the beginning and end of the data sequence, however, and these artifacts are due to the moving window problem noted above for the first K and the last K elements of the original data sequence. In particular, the filter implementation in the savgol procedure does not apply the sequence extension procedure discussed above, and this fact is responsible for these artifacts appearing at the beginning and end of the smoothed data sequence.

It is extremely easy to correct this problem, adopting the same philosophy the package uses for the outlierMAD procedure: simply retain the first and last K elements of the original sequence unmodified. The procedure SGwrapper listed below does this after the fact, calling the savgol procedure and then replacing the first and last K elements of the filtered sequence with the original sequence values:

SGwrapper <- function(x,K,forder=4,dorder=0){

n = length(x)

fl = 2*K+1

y = savgol(x,fl,forder,dorder)

if (dorder == 0){

y[1:K] = x[1:K]

y[(n-K):n] = x[(n-K):n]

}

else{

y[1:K] = 0

y[(n-K):n] = 0

}

Before showing the results obtained with this procedure, it is important to note two points. First, the moving window width parameter fl required for the savgol procedure corresponds to fl = 2K+1 for a half-width parameter K. The procedure SGwrapper instead requires K as its passing parameter, constructing fl from this value of K. Second, note that in addition to serving as a smoother, the Savitzky-Golay filter family can also be used to estimate derivatives (this is tricky since differentiation filters are incredible noise amplifiers, but I’ll talk more about that in another post). In the savgol procedure, this is accomplished by specifying the parameter dorder, which has a default value of zero (implying smoothing), but which can be set to 1 to estimate the first derivative of a sequence, 2 for the second derivative, etc. In these cases, replacing the first and last K elements of the filtered sequence with the original data sequence elements is not reasonable: in the absence of any other knowledge, a better default derivative estimate is zero, and the SGwrapper procedure listed above does this.

The four plots shown above illustrate the differences between the original savgol procedure (the left-hand plots) and those obtained with the SGwrapper procedure listed above (the right-hand plots). In all cases, the data sequence used to generate these plots was the physical property data sequence cleaned using the Hampel filter with t0 = 3. The upper left plot repeats the lower of the two previous plots, corresponding to the savgol smoother output, while the upper right plot applies the SGwrapper function to remove the artifacts at the beginning and end of the smoothed data sequence. Similarly, the lower two plots give the corresponding second-derivative estimates, obtained by applying the savgol procedure with fl = 11 and dorder = 2 (lower left plot) or the SGwrapper procedure with K = 5 and dorder = 2 (lower right plot).

Cleaning time-series and other data streams

2011-11-27T08:37:00.000-08:00

The need to analyze time-series or other forms of streaming data arises frequently in many different application areas. Examples include economic time-series like stock prices, exchange rates, or unemployment figures, biomedical data sequences like electrocardiograms or electroencephalograms, or industrial process operating data sequences like temperatures, pressures or concentrations. As a specific example, the figure below shows four data sequences: the upper two plots represent hourly physical property measurements, one made at the inlet of a product storage tank (the left-hand plot) and the other made at the same time at the outlet of the tank (the right-hand plot). The lower two plots in this figure show the results of applying the data cleaning filter outlierMAD from the R package pracma discussed further below. The two main points of this post are first, that isolated spikes like those seen in the upper two plots at hour 291 can badly distort the results of an otherwise reasonable time-series characterization, and second, that the simple moving window data cleaning filter described here is often very effective in removing these artifacts.

This example is discussed in more detail in Section 8.1.2 of my book Discrete-Time Dynamic Models, but the key observations here are the following. First, the large spikes seen in both of the original data sequences were caused by the simultaneous, temporary loss of both measurements and the subsequent coding of these missing values as zero by the data collection system. The practical question of interest was to determine how long, on average, the viscous, polymeric material being fed into and out of the product storage tank was spending there. A standard method for addressing such questions is the use of cross-correlation analysis, where the expected result is a broad peak like the heavy dashed line in the plot shown below. The location of this peak provides an estimate of the average time spent in the tank, which is approximately 21 hours in this case, as indicated in the plot. This result was about what was expected, and it was obtained by applying standard cross-correlation analysis to the cleaned data sequences shown in the bottom two plots above. The lighter solid curve in the plot below shows the results of applying exactly the same analysis, but to the original data sequences instead of the cleaned data sequences. This dramatically different plot suggests that the material is spending very little time in the storage tank: accepted uncritically, this result would imply severe fouling of the tank, suggesting a need to shut the process down and clean out the tank, an expensive and labor-intensive proposition. The main point of this example is that the difference in these two plots is entirely due to the extreme data anomalies present in the original time-series. Additional examples of problems caused by time-series outliers are discussed in Section 4.3 of my book Mining Imperfect Data.

One of the primary features of the analysis of time-series and other streaming data sequences is the need for local data characterizations. This point is illustrated in the plot below, which shows the first 200 observations of the storage tank inlet data sequence discussed above. All of these observations but one are represented as open circles in this plot, but the data point at k = 110 is shown as a solid circle, to emphasize how far it lies from its immediate neighbors in the data sequence. It is important to note that this point is not anomalous with respect to the overall range of this data sequence – it is, for example, well within the normal range of variation seen for the points from about k = 150 to k = 200 – but it is clearly anomalous with respect to those points that immediately precede and follow it. A general strategy for automatically detecting and removing such spikes from a data sequence like this one is to apply a moving window data cleaning filter which characterizes each data point with respect to a local neighborhood of prior and subsequent samples. That is, for each data point k in the original data sequence, this type of filter forms a cleaned data estimate based on some number J of prior data values (i.e., points k-J through k-1 in the sequence) and, in the simplest implementations, the same number of subsequent data values (i.e., points k+1 through k+J in the sequence).

The specific data cleaning filter considered here is the Hampel filter, which applies the Hampel identifier discussed in Chapter 7 of Exploring Data in Engineering, the Sciences and Medicine to this moving data window. If the k^th data point is declared to be an outlier, it is replaced by the median value computed from this data window; otherwise, the data point is not modified. The results of applying the Hampel filter with a window width of J = 5 to the above data sequence are shown in the plot below. The effect is to modify three of the original data points – those at k = 43, 110, and 120 – and the original values of these modified points are shown as solid circles at the appropriate locations in this plot. It is clear that the most pronounced effect of the Hampel filter is to remove the local outlier indicated in the above figure and replace it with a value that is much more representative of the other data points in the immediate vicinity.

As I noted above, the Hampel filter implementation used here is that available in the R package pracma as procedure outlierMAD. I will discuss this R package in more detail in my next post, but for those seeking a more detailed discussion of the Hampel filter in the meantime, one is freely available on-line in the form of an EDN article I wrote in 2002, Scrub data with scale-invariant nonlinear digital filters. Also, comparisons with alternatives like the standard median filter (generally too aggressive, introducing unwanted distortion into the “cleaned” data sequence) and the center-weighted median filter (sometimes quite effective) are presented in Section 4.2 of the book Mining Imperfect Data mentioned above.

Harmonic means, reciprocals, and ratios of random variables

2011-11-11T14:16:00.000-08:00

In my last few posts, I have considered “long-tailed” distributions whose probability density decays much more slowly than standard distributions like the Gaussian. For these slowly-decaying distributions, the harmonic mean often turns out to be a much better (i.e., less variable) characterization than the arithmetic mean, which is generally not even well-defined theoretically for these distributions. Since the harmonic mean is defined as the reciprocal of the mean of the reciprocal values, it is intimately related to the reciprocal transformation. The main point of this post is to show how profoundly the reciprocal transformation can alter the character of a distribution, for better or worse. One way that reciprocal transformations sneak into analysis results is through attempts to characterize ratios of random numbers. The key issue underlying all of these ideas is the question of when the denominator variable in either a reciprocal transformation or a ratio exhibits non-negligible probability in a finite neighborhood of zero. I discuss transformations in Chapter 12 of Exploring Data in Engineering, the Sciences and Medicine, with a section (12.7) devoted to reciprocal transformations, showing what happens when we apply them to six different distributions: Gaussian, Laplace, Cauchy, beta, Pareto, and lognormal.

In the general case, if a random variable x has the density p(x), the distribution g(y) of the reciprocal y = 1/x has the density:

g(y) = p(1/y)/y²

As I discuss in greater detail in Exploring Data, the consequence of this transformation is typically (though not always) to convert a well-behaved distribution into a very poorly behaved one. As a specific example, the plot below shows the effect of the reciprocal transformation on a Gaussian random variable with mean 1 and standard deviation 2. The most obvious characteristic of this transformed distribution is its strongly asymmetric, bimodal character, but another non-obvious consequence of the reciprocal transformation is that it takes a distribution that is completely characterized by its first two moments into a new distribution with Cauchy-like tails, for which none of the integer moments exist.

The implications of the reciprocal transformation for many other distributions are equally non-obvious. For example, both the badly-behaved Cauchy distribution (no moments exist) and the well-behaved lognormal distribution (all moments exist, but interestingly, do not completely characterize the distribution, as I have discussed in a previous post) are invariant under the reciprocal transformation. Also, applying the reciprocal transformation to the long-tailed Pareto type I distribution (which exhibits few or no finite moments, depending on its tail decay rate) yields a beta distribution, all of whose moments are finite. Finally, it is worth noting that the invariance of the Cauchy distribution under the reciprocal transformation lies at the heart of the following result, presented in the book Continuous Univariate Distributions by Johnson, Kotz, and Balakrishnan (Volume 1, 2^nd edition, Wiley, 1994, page 319). They note that if the density of x is positive, continuous, and differentiable at x = 0 – all true for the Gaussian case – the distribution of the harmonic mean of N samples approaches a Cauchy limit as N becomes infinitely large.

As noted above, the key issue responsible for the pathological behavior of the reciprocal transformation is the question of whether the original data distribution exhibits nonzero probability of taking on values within a neighborhood around zero. In particular, note that if x can only assume values larger than some positive lower limit L, it follows that 1/x necessarily lies between 0 and 1/L, which is enough to guarantee that all moments of the transformed distribution exist. For the Gaussian distribution, even if the mean is large enough and the standard deviation is small enough that the probability of observing values less than some limit L > 0 is negligible, the fact that this probability is not zero means that the moments of any reciprocally-transformed Gaussian distribution are not finite. As a practical matter, however, reciprocal transformations and related characterizations – like harmonic means and ratios – do become better-behaved as the probability of observing values near zero become negligibly small.

To see this point, consider two reciprocally-transformed Gaussian examples. The first is the one considered above: the reciprocal transformation of a Gaussian random variable with mean 1 and standard deviation 2. In this case, the probability that x assumes values smaller than or equal to zero is non-negligible. Specifically, this probability is simply the cumulative distribution function for the distribution evaluated at zero, easily computed in R as approximately 31%:

> pnorm(0,mean=1,sd=2)

[1] 0.3085375

In contrast, for a Gaussian random variable with mean 1 and standard deviation 0.1, the corresponding probability is negligibly small:

> pnorm(0,mean=1,sd=0.1)

[1] 7.619853e-24

If we consider the harmonic means of these two examples, we see that the first one is horribly behaved, as all of the results presented here would lead us to expect. In fact, the qqPlot command in the car package in R allows us to compute quantile-quantile plots for the Student’s t-distribution with one degree of freedom, corresponding to the Cauchy distribution, yielding the plot shown below. The Cauchy-like tail behavior expected from the results presented by Johnson, Kotz and Balakrishnan is seen clearly in this Cauchy Q-Q plot, constructed from 1000 harmonic means, each computed from statistically independent samples drawn from a Gaussian distribution with mean 1 and standard deviation 2. The fact that almost all of the observations fall within the – very wide – 95% confidence interval around the reference line suggest that the Cauchy tail behavior is appropriate here.

To further confirm this point, compare the corresponding normal Q-Q plot for the same sequence of harmonic means, shown below. There, the extreme non-Gaussian character of these harmonic means is readily apparent from the pronounced outliers evident in both the upper and lower tails.

In marked contrast, for the second example with the mean of 1 as before but the much smaller standard deviation of 0.1, the harmonic mean is much better behaved, as the normal Q-Q plot below illustrates. Specifically, this plot is identical in construction to the one above, except it was computed from samples drawn from the second data distribution. Here, most of the computed harmonic mean values fall within the 95% confidence limits around the Gaussian reference line, suggesting that it is not unreasonable in practice to regard these values as approximately normally distributed, in spite of the pathologies of the reciprocal transformation.

One reason the reciprocal transformation is important in practice – particularly in connection with the Gaussian distribution – is that the desire to characterize ratios of uncertain quantities does arise from time to time. In particular, if we are interested in characterizing the ratio of two averages, the Central Limit Theorem would lead us to expect that, at least approximately, this ratio should behave like the ratio of two Gaussian random variables. If these component averages are statistically independent, the expected value of the ratio can be re-written as the product of the expected value of the numerator average and the expected value of the reciprocal of the denominator average, leading us directly to the reciprocal Gaussian transformation discussed here. In fact, if these two averages are both zero mean, it is a standard result that the ratio has a Cauchy distribution (this result is presented in the same discussion from Johnson, Kotz and Balakrishnan noted above). As in the second harmonic mean example presented above, however, it turns out to be true that if the mean and standard deviation of the denominator variable are such that the probability of a zero or negative denominator are negligible, the distribution of the ratio may be approximated reasonably well as Gaussian. A very readable and detailed discussion of this fact is given in the paper by George Marsaglia in the May 2006 issue of Journal of Statistical Software.

Finally, it is important to note that the “reciprocally-transformed Gaussian distribution” I have been discussing here is not the same as the inverse Gaussian distribution, to which Johnson, Kotz and Balakrishnan devote a 39-page chapter (Chapter 15). That distribution takes only positive values and exhibits moments of all orders, both positive and negative, and as a consequence, it has the interesting characteristic that it remains well-behaved under reciprocal transformations, in marked contrast to the Gaussian case.

The Zipf and Zipf-Mandelbrot distributions

2011-10-23T13:31:00.000-07:00

In my last few posts, I have been discussing some of the consequences of the slow decay rate of the tail of the Pareto type I distribution, along with some other, closely related notions, all in the context of continuously distributed data. Today’s post considers the Zipf distribution for discrete data, which has come to be extremely popular as a model for phenomena like word frequencies, city sizes, or sales rank data, where the values of these quantities associated with randomly selected samples can vary by many orders of magnitude.

More specifically, the Zipf distribution is defined by a probability p_i of observing the i^th element of an infinite sequence of objects in a single random draw from that sequence, where the probability is given by:

p_i = A/i^a

Here, a is a positive number greater than 1 that determines the rate of the distribution’s tail decay, and A is a normalization constant, chosen so that these probabilities sum to 1. Like the continuous-valued Pareto type I distribution, the Zipf distribution exhibits a “long tail,” meaning that its tail decays slowly enough that in a random sample of objects O_i drawn from a Zipf distribution, some very large values of the index i will be observed, particularly for relatively small values of the exponent a. In one of the earliest and most common applications of the Zipf distribution, the objects considered represent words in a document and i represents their rank, ranging from most frequent (for i = 1) to rare (for large i ). In a more business-oriented application, the objects might be products for sale (e.g., books listed on Amazon), with the index i corresponding to their sales rank. For a fairly extensive collection of references to many different applications of the Zipf distribution, the website (originally) from Rockefeller University is an excellent source.

In Exploring Data in Engineering, the Sciences, and Medicine, I give a brief discussion of both the Zipf distribution and the closely related Zipf-Mandelbrot distribution discussed by Beniot Mandelbrot in his book The Fractal Geometry of Nature. The probabilities defining this distribution may be parameterized in several ways, and the one given in Exploring Data is:

p_i = A/(1+Bi)^a

where again a is an exponent that determines the rate at which the tail of the distribution decays, and B is a second parameter with a value that is strictly positive but no greater than 1. For both the Zipf distribution and the Zipf-Mandelbrot distribution, the exponent a must be greater than 1 for the distribution to be well-defined, it must be greater than 2 for the mean to be finite, and it must be greater than 3 for the variance to be finite.

So far, I have been unable to find an R package that supports the generation of random samples drawn from the Zipf distribution, but the package zipfR includes the command rlnre, which generates random samples drawn from the Zipf-Mandelbrot distribution. As I noted, this distribution can be parameterized in several different ways and, as Murphy’s law would have it, the zipfR parameterization is not the same as the one presented above and discussed in Exploring Data. Fortunately, the conversion between these parameters is simple. The zipfR package defines the distribution in terms of a parameter alpha that must lie strictly between 0 and 1, and a second parameter B that I will call B_zipfR to avoid confusion with the parameter B in the above definition. These parameters are related by:

alpha = 1/a and B_zipfR = (a-1) B

Since the a parameter (and thus the alpha parameter in the zipfR package) determines the tail decay rate of the distribution, it is of the most interest here, and the rest of this post will focus on three examples: a = 1.5 (alpha = 2/3), for which both the distribution’s mean and variance are infinite, a = 2.5 (alpha = 2/5), for which the mean is finite but the variance is not, and a = 3.5 (alpha = 2/7), for which both the mean and variance are finite. The value of the parameter B in the Exploring Data definition of the distribution will be fixed at 0.2 in all of these examples, corresponding to values of B_zipfR = 0.1, 0.3, and 0.5 for the three examples considered here.

To generate Zipf-Mandelbrot random samples, the zipfR package uses the procedure rlnre in conjunction with the procedure lnre (the abbreviation “lnre” stands for “large number of rare events” and it represents a class of data models that includes the Zipf-Mandelbrot distribution). Specifically, to generate a random sample of size N = 100 for the first case considered here, the following R code is executed:

> library(zipfR)
> ZM = lnre(“zm”, alpha = 2/3, B = 0.1)
> zmsample = rlnre(ZM, n=100)

The first line loads the zipfR library (which must first be installed, of course, using the install.packages command), the second line invokes the lnre command to set up the distribution with the desired parameters, and the last line invokes the rlnre command to generate 100 random samples from this distribution. (As with all R random number generators, the set.seed command should be used first to initialize the random number generator seed if you want to get repeatable results; for the results presented here, I used set.seed(101).) The sample returned by the rlnre command is a vector of 100 observations, which have the “factor” data type, although their designations are numeric (think of the factor value “1339” as meaning “1 sample of object number 1339”). In the results I present here, I have converted these factor responses to numerical ones so I can interpret them as numerical ranks. This conversion is a little subtle: simply converting from factor to numeric values via something like “zmnumeric = as.numeric(zmsample)” almost certainly doesn’t give you what you want: this will convert the first-ocurring factor value (which has a numeric label, say “1339”) into the number 1, convert the second-occurring value (since this is a random sequence, this might be “73”) into the number 2, etc. To get what you want (e.g., the labels “1339” and “73” assigned to the numbers 1339 and 73, respectively), you need to first convert the factors in zmsample into characters and then convert these characters into numeric values:

zmnumeric = as.numeric(as.character(zmsample))

The three plots below show random samples drawn from each of the three Zipf-Mandelbrot distributions considered here. In all cases, the y-axis corresponds to the number of times the object labeled i was observed in a random sample of size N = 100 drawn from the distribution with the indicated exponent. Since the range of these indices can be quite large in the slowly-decaying members of the Zipf-Mandelbrot distribution family, the plots are drawn with logarithmic x-axes, and to facilitate comparisons, the x-axes have the same range in all three plots, as do the y-axes. In all three plots, object i = 1 occurs most often – about a dozen times in the top plot, two dozen times in the middle plot, and three dozen times in the bottom plot – and those objects with larger indices occur less frequently. The major difference between these three examples lies in the largest indices of the objects seen in the samples: we never see an object with index greater than 50 in the bottom plot, we see only two such objects in the middle plot, while more than a third of the objects in the top plot meet this condition, with the most extreme object having index i = 115,116.

As in the case of the Pareto type I distributions I discussed in several previous posts – which may be regarded as the continuous analog of the Zipf distribution – the mean is generally not a useful characterization for the Zipf distribution. This point is illustrated in the boxplot comparison presented below, which summarizes the means computed from 1000 statistically independent random samples drawn from each of the three distributions considered here, where the object labels have been converted to numerical values as described above. Thus, the three boxplots on the left represent the means – note the logarithmic scale on the y-axis – of these index values i generated for each random sample. The extreme variability seen for Case 1 (a = 1.5) reflects the fact that neither the mean nor the variance are finite for this case, and the consistent reduction in the range of variability for Cases 2 (a = 2.5, finite mean but infinite variance) and 3 (a = 3.5, finite mean and variance) reflects the “shortening tail” of this distribution with increasing exponent a. As I discussed in my last post, a better characterization than the mean for distributions like this is the “95% tail length,” corresponding to the 95% sample quantile. Boxplots summarizing these values for the three distributions considered here are shown to the right of the dashed vertical line in the plot below. In each case, the range of variation seen here is much less extreme for the 95% tail length than it is for the mean, supporting the idea that this is a better characterization for data described by Zipf-like discrete distributions.

Other alternatives to the (arithmetic) mean that I discussed in conjunction with the Pareto type I distribution were the sample median, the geometric mean, and the harmonic mean. The plot below compares these four characterizations for 1000 random samples, each of size N = 100, drawn from the Zipf-Mandelbrot distribution with a = 3.5 (the third case), for which the mean is well-defined. Even here, it is clear that the mean is considerably more variable than these other three alternatives.

Finally, the plot below shows boxplot comparisons of these alternative characterizations – the median, the geometric mean, and the harmonic mean – for all three of the distributions considered here. Not surprisingly, Case 1 (a = 1.5) exhibits the largest variability seen for all three characterizations, but the harmonic mean is much more consistent for this case than either the geometric mean or the median. In fact, the same observation holds – although less dramatically – for Case 2 (a = 2.5), and the harmonic mean appears more consistent than the geometric mean for all three cases. This observation is particularly interesting in view of the connection between the harmonic mean and the reciprocal transformation, which I will discuss in more detail next time.

Is the “Long Tail” a Useless Concept?

2011-09-28T15:11:00.000-07:00

In response to my last post, “The Long Tail of the Pareto Distribution,” Neil Gunther had the following comment:

“Unfortunately, you've fallen into the trap of using the ‘long tail’ misnomer. If you think about it, it can't possibly be the length of the tail that sets distributions like Pareto and Zipf apart; even the negative exponential and Gaussian have infinitely long tails.”

He goes on to say that the relevant concept is the “width” or the “weight” of the tails that is important, and that a more appropriate characterization of these “Long Tails” would be “heavy-tailed” or “power-law” distributions.

Neil’s comment raises an important point: while the term “long tail” appears a lot in both the on-line and hard-copy literature, it is often somewhat ambiguously defined. For example, in his book, The Long Tail, Chris Anderson offers the following description (page 10):

“In statistics, curves like that are called ‘long-tailed distributions’ because the tail of the curve is very long relative to the head.”

The difficulty with this description is that it is somewhat ambiguous since it says nothing about how to measure “tail length,” forcing us to adopt our own definitions. It is clear from Neil’s comments that the definition he adopts for “tail length” is the width of the distribution’s support set. Under this definition, the notion of a “long-tailed distribution” is of extremely limited utility: the situation is exactly as Neil describes it, with “long-tailed distributions” corresponding to any distribution with unbounded support, including both distributions like the Gaussian and gamma distribution where the mean is a reasonable characterization, and those like the Cauchy and Pareto distribution where the mean doesn’t even exist.

The situation is analogous to that of confidence intervals, which characterize the uncertainty inherited by any characterization computed from a collection of uncertain (i.e., random) data values. As a specific example consider the mean: the sample mean is the arithmetic average of N observed data samples, and it is generally intended as an estimate of the population mean, defined as the first moment of the data distribution. A q% confidence interval around the sample mean is an interval that contains the population mean with probability at least q%. These intervals can be computed in various ways for different data characterizations, but the key point here is that they are widely used in practice, with the most popular choices being the 90%, 95% and 99% confidence intervals, which necessarily become wider as this percentage q increases. (For a more detailed discussion of confidence intervals, refer to Chapter 9 of Exploring Data in Engineering, the Sciences, and Medicine.) We can, in principle, construct 100% confidence intervals, but this leads us directly back to Neil’s objection: the 100% confidence interval for the mean is entire support set of the distribution (e.g., for the Gaussian distribution, this 100% confidence interval is the whole real line, while for any gamma distribution, it is the set of all positive numbers). These observations suggest the following notion of “tail length” that addresses Neil’s concern while retaining the essential idea of interest in the business literature: we can compare the “q% tail length” of different distributions for some q less than 100.

In particular, consider the case of J-shaped distributions, defined as those like the Pareto type I distribution whose distribution p(x) decays monotonically with increasing x, approaching zero as x goes to infinity. The plot below shows two specific examples to illustrate the idea: the solid line corresponds to the (shifted) exponential distribution:

p(x) = e^–(x-1)

for all x greater than or equal to 1 and zero otherwise, while the dotted line represents the Pareto type I distribution with location parameter k = 1 and shape parameter a = 0.5 discussed in my last post. Initially, as x increases from 1, the exponential density is greater than the Pareto density, but for x larger than about 3.5, the opposite is true: the exponential distribution rapidly becomes much smaller, reflecting its much more rapid rate of tail decay.

For these distributions, define the q% tail length to be the distance from the minimum possible value of x (the “head” of the distribution; here, x = 1) to the point in the tail where the cumulative probability reaches q% (i.e., the value x_q where x < x_q with probability q%). In practical terms, the q% tail length tells us how far out we have to go in the tail to account for q% of the possible cases. In R, this value is easy to compute using the quantile function included in most families of available distribution functions. As a specific example, for the Pareto type I distribution, the function qparetoI in the VGAM package gives us the desired quantiles for the distribution with specified values of the parameters k (designated “scale” in the qparetoI call) and a (designated “shape” in the qparetoI call). Thus, for the case k = 1 and a = 0.5 (i.e., the dashed curve in the above plot), the “90% tail length” is given by:

> qparetoI(p=0.9,scale=1,shape=0.5)

[1] 100

For comparison, the corresponding shifted exponential distribution has the 90% tail length given by:

> 1 + qexp(p = 0.9)

[1] 3.302585

(Note that here, I added 1 to the exponential quantile to account for the shift in its domain from “all positive numbers” – the domain for the standard exponential distribution – to the shifted domain “all numbers greater than 1”.) Since these 90% tail lengths differ by a factor of 30, they provide a sound basis for declaring the Pareto type I distribution to be “longer tailed” than the exponential distribution.

These results also provide a useful basis for assessing the influence of the decay parameter a for the Pareto distribution. As I noted last time, two of the examples I considered did not have finite means (a = 0.5 and 1.0), and none of the four had finite variances (i.e., also a = 1.5 and 2.0), rendering moment characterizations like the mean and standard deviation fundamentally useless. Comparing the 90% tail lengths for these distributions, however, leads to the following results:

a = 0.5: 90% tail length = 100.000

a = 1.0: 90% tail length = 10.000

a = 1.5: 90% tail length = 4.642

a = 2.0: 90% tail length = 3.162

It is clear from these results that the shape parameter a has a dramatic effect on the 90% tail length (in fact, on the q% tail length for any q less than 100). Further, note that the 90% tail length for the Pareto type I distribution with a = 2.0 is actually a little bit shorter than that for the exponential distribution. If we move further out into the tail, however, this situation changes. As a specific example, suppose we compare the 98% tail lengths. For the exponential distribution, this yields the value 4.912, while for the four Pareto shape parameters we have the following results:

a = 0.5: 98% tail length = 2,500.000

a = 1.0: 98% tail length = 50.000

a = 1.5: 98% tail length = 13.572

a = 2.0: 98% tail length = 7.071

This value (i.e., the 98% tail length) seems a particularly appropriate choice to include here since in his book, The Long Tail, Chris Anderson notes that his original presentations on the topic were entitled “The 98% Rule,” reflecting the fact that he was explicitly considering how far out you had to go into the tail of a distribution of goods (e.g., the books for sale by Amazon) to account for 98% of the sales.

Since this discussion originally began with the question, “when are averages useless?” it is appropriate to note that, in contrast to the much better-known average, the “q% tail length” considered here is well-defined for any proper distribution. As the examples discussed here demonstrate, this characterization also provides a useful basis for quantifying the “Long Tail” behavior that is of increasing interest in business applications like Internet marketing. Thus, if we adopt this measure for any q value less than 100%, the answer to the title question of this post is, “No: The Long Tail is a useful concept.”

The downside of this minor change is that – as the results shown here illustrate – the results obtained using the q% tail length depend on the value of q we choose. In my next post, I will explore the computational issues associated with that choice.

The Long Tail of the Pareto Distribution

2011-09-17T09:54:00.000-07:00

In my last two posts, I have discussed cases where the mean is of little or no use as a data characterization. One of the specific examples I discussed last time was the case of the Pareto type I distribution, for which the density is given by:

p(x) = ak^a/x^a+1

defined for all x > k, where k and a are numeric parameters that define the distribution. In the example I discussed last time, I considered the case where a = 1.5, which exhibits a finite mean (specifically, the mean is 3 for this case), but an infinite variance. As the results I presented last time demonstrated, the extreme data variability of this distribution renders the computed mean too variable to be useful. Another reason this distribution is particularly interesting is that it exhibits essentially the same tail behavior as the discrete Zipf distribution; there, the probability that a discrete random variable x takes its i^th value is:

p_i = A/i^c,

where A is a normalization constant and c is a parameter that determines how slowly the tail decays. This distribution was originally proposed to characterize the frequency of words in long documents (the Zipf-Estoup law), it was investigated further by Zipf in the mid-twentieth century in a wide range of applications (e.g., the distributions of city sizes), and it has become the subject of considerable recent attention as a model for “long-tailed” business phenomena (for a non-technical introduction to some of these business phenomena, see the book by Chris Anderson, The Long Tail). I will discuss the Zipf distribution further in a later post, but one of the reasons for discussing the Pareto type I distribution first is that since it is a continuous distribution, the math is easier, meaning that more characterization results are available for the Pareto distribution.

The mean of the Pareto type I distribution is:

Mean = ak/(a-1),

provided a > 1, and the variance of the distribution is finite only if a > 2. Plots of the probability density defined above for this distribution are shown above, for k = 1 in all cases, and with a taking the values 0.5, 1.0, 1.5, and 2.0. (This is essentially the same plot as Figure 4.17 in Exploring Data in Engineering, the Sciences, and Medicine, where I give a brief description of the Pareto type I distribution.) Note that all of the cases considered here are characterized by infinite variance, while the first two (a = 0.5 and 1.0) are also characterized by infinite means. As the results presented below emphasize, the mean represents a very poor characterization in practice for data drawn from any of these distributions, but there are alternatives, including the familiar median that I have discussed previously, along with two others that are more specific to the Pareto type I distribution: the geometric mean and the harmonic mean.

The plot below emphasizes the point made above about the extremely limited utility of the mean as a characterization of Pareto type I data, even in cases where it is theoretically well-defined. Specifically, this plot compares the four characterizations I discuss here – the mean (more precisely known as the “arithmetic mean” to distinguish it from the other means considered here), the median, the geometric mean, and the harmonic mean – for 1000 statistically independent Pareto type I data sequences, each of length N = 400, with parameters k = 1 and a = 2.0. For this example, the mean is well-defined (specifically, it is equal to 2), but compared with the other data characterizations, its variability is much greater, reflecting the more serious impact of this distribution’s infinite variance on the mean than on these other data characterizations.

To give a more complete view of the extreme variability of the arithmetic mean, boxplots of 1000 statistically independent samples drawn from all four of the Pareto type I distribution examples plotted above are shown in the boxplots below. As before, each sample is of size N = 400 and the parameter k has the value 1, but here the computed arithmetic means are shown for the parameter values a = 0.5, 1.0, 1.5, and 2.0; note the log scale used here because the range of computed means is so large. For the first two of these examples, the population mean does not exist, so it is not surprising that the computed values span such an enormous range, but even when the mean is well-defined, the influence of the infinite variance of these cases is clearly evident. It may be argued that infinite variance is an extreme phenomenon, but it is worth emphasizing here that for the specific “long tail” distributions popular in many applications, the decay rate is sufficiently slow for the variance – and sometimes even the mean – to be infinite, as in these examples.

As I have noted several times in previous posts, the median is much better behaved than the mean, so much so that it is well-defined for any proper distribution. One of the advantages of the Pareto type I distribution is that the form of the density function is simple enough that the median of the distribution can be computed explicitly from the distribution parameters. This result is given in the fabulous book by Johnson, Kotz and Balakrishnan that I have mentioned previously, which devotes an entire chapter (Chapter 20) to the Pareto family of distributions. Specifically, the median of the Pareto type I distribution with parameters k and a is given by:

Median = 2^1/ak

Thus, for the four examples considered here, the median values are 4.0 (for a = 0.5), 2.0 (for a = 1.0), 1.587 (for a = 1.5), and 1.414 (for a = 2.0). Boxplot summaries for the same 1000 random samples considered above are shown in the plot below, which also includes horizontal dotted lines at these theoretical median values for the four distributions. The fact that these lines correspond closely with the median lines in the boxplots gives an indication that the computed median is, on average, in good agreement with the correct values it is attempting to estimate. As in the case of the arithmetic means, the variability of these estimates decreases monotonically as a increases, corresponding to the fact that the distribution becomes generally better-behaved as the a parameter increases.

The geometric mean is an alternative characterization to the more familiar arithmetic mean, one that is well-defined for any sequence of positive numbers. Specifically, the geometric mean of N positive numbers is defined as the N^th root of their product. Equivalently, the geometric mean may be computed by exponentiating the arithmetic average of the log-transformed values. In the case of the Pareto type I distribution, the utility of the geometric mean is closely related to the fact that the log transformation converts a Pareto-distributed random variable into an exponentially distributed one, a point that I will discuss further in a later post on data transformations. (These transformations are the topic of Chapter 12 of Exploring Data, where I briefly discuss both the logarithmic transformation on which the geometric mean is based and the reciprocal transformation on which the harmonic mean is based, described next.) The key point here is that the following simple expression is available for the geometric mean of the Pareto type I distribution (Johnson, Kotz, and Balakrishnan, page 577):

Geometric Mean = k exp(1/a)

For the four specific examples considered here, these geometric mean values are approximately 7.389 (for a = 0.5), 2.718 (for a = 1.0), 1.948 (for a = 1.5), and 1.649 (for a = 2.0). The boxplots shown below summarize the range of variation seen in the computed geometric means for the same 1000 statistically independent samples considered above. Again, the horizontal dotted lines indicate the correct values for each distribution, and it may be seen that the computed values are in good agreement, on average. As before, the variability of these computed values decreases with increasing a values as the distribution becomes better-behaved.

The fourth characterization considered here is the harmonic mean, again appropriate to positive values, and defined as the reciprocal of the average of the reciprocal data values. In the case of the geometric mean just discussed, the log transformation on which it is based is often useful in improving the distributional character of data values that span a wide range. In the case of the Pareto type I distribution – and a number of others – the reciprocal transformation on which the harmonic mean is based also improves the behavior of the data distribution, but this is often not the case. In particular, reciprocal transformations often make the character of a data distribution much worse: applied to the extremely well-behaved standard uniform distribution, it yields the Pareto type I distribution with a = 1, for which none of the integer moments exist; similarly, applied to the Gaussian distribution, the reciprocal transformation yields a result that is both infinite variance and bimodal. (A little thought suggests that the reciprocal transformation is inappropriate for the Gaussian distribution because it is not strictly positive, but normality is a favorite working assumption, sometimes applied to the denominators of ratios, leading to a number of theoretical difficulties. I will have more to say about that in a future post.) For the case of the Pareto type I distribution, the reciprocal transformation converts it into the extremely well-behaved beta distribution, and the harmonic mean has the following simple expression:

Harmonic mean = k(1 + a^-1)

For the four examples considered here, this expression yields harmonic mean values of 3 (for a = 0.5), 2 (for a = 1.0), 1.667 (for a = 1.5), and 1.5 (for a = 2.0). Boxplot summaries of the computed harmonic means for the 1000 simulations of each case considered previously are shown below, again with dotted horizontal lines at the theoretical values for each case. As with both the median and the geometric mean, it is clear from these plots that the computed values are correct on average, and their variability decreases with increasing values of the a parameter.

The key point of this post has been to show that, while averages are not suitable characterizations for “long tailed” phenomena that are becoming an increasing subject of interest in many different fields, useful alternatives do exist. For the case of the Pareto type I distribution considered here, these alternatives include the popular median, along with the somewhat less well-known geometric and harmonic means. In an upcoming post, I will examine the utility of these characterizations for the Zipf distribution.

Some Additional Thoughts on Useless Averages

2011-08-27T13:46:00.000-07:00

In my last post, I described three situations where the average of a sequence of numbers is not representative enough to be useful: in the presence of severe outliers, in the face of multimodal data distributions, and in the face of infinite-variance distributions. The post generated three interesting comments that I want to respond to here.

First and foremost, I want to say thanks to all of you for giving me something to think about further, leading me in some interesting new directions. First, chrisbeeleyimh had the following to say:

“I seem to have rather abandoned means and medians in favor of drawing the distribution all the time, which baffles my colleagues somewhat.”

Chris also maintains a collection of data examples where the mean is the same but the shape is very different. In fact, one of the points I illustrate in Section 4.4.1 of Exploring Data in Engineering, the Sciences, and Medicine is that there are cases where not only the means but all of the moments (i.e., variance, skewness, kurtosis, etc.) are identical but the distributions are profoundly different. A specific example is taken from the book Counterexamples in Probability, 2nd Edition by J.M. Stoyanov, who shows that if the lognormal density is multiplied by the following function:

f(x) = 1 + A sin(2 pi ln x),

for any constant A between -1 and +1, the moments are unchanged. The character of the distribution is changed profoundly, however, as the following plot illustrates (this plot is similar to Fig. 4.8 in Exploring Data, which shows the same two distributions, but for A = 0.5 instead of A = 0.9, as shown here). To be sure, this behavior is pathological – distributions that have finite support, for example, are defined uniquely by their complete set of moments – but it does make the point that moment characterizations are not always complete, even if an infinite number of them are available. Within well-behaved families of distributions (such as the one proposed by Karl Pearson in 1895), a complete characterization is possible on the basis of the first few moments, which is one reason for the historical popularity of the method of moments for fitting data to distributions. It is important to recognize, however, that moments do have their limitations and that the first moment alone – i.e., the mean by itself – is almost never a complete characterization. (I am forced to say “almost” here because if we impose certain very strong distributional assumptions – e.g., the Poisson or binomial distributions – the specific distribution considered may be fully characterized by its mean. This begs the question, however, of whether this distributional assumption is adequate. My experience has been that, no matter how firmly held the belief in a particular distribution is, exceptions do arise in practice … overdispersion, anyone?)

The plot below provides a further illustration of the inadequacy of the mean as a sole data characterization, comparing four different members of the family of beta distributions. These distributions – in the standard form assumed here – describe variables whose values range from 0 to 1, and they are defined by two parameters, p and q, that determine the shape of the density function and all moments of the distribution. The mean of the beta distribution is equal to p/(p+q), so if p = q – corresponding to the class of symmetric beta distributions – the mean is ½, regardless of the common value of these parameters. The four plots below show the corresponding distributions when both parameters are equal to 0.5 (upper left, the arcsin distribution I discussed last time), 1.0 (upper right, the uniform distribution), 1.5 (lower left), and 8.0 (lower right).

The second comment on my last post was from Efrique, who suggested the Student’s t-distribution with 2 degrees of freedom as a better infinite-variance example than the Cauchy example I used (corresponding to Student’s t-distribution with one degree of freedom), because the first moment doesn’t even exist for the Cauchy distribution (“there’s nothing to converge to”). The figure below expands the boxplot comparison I presented last time, comparing the means, medians, and modes (from the modeest package), for both of these infinite-variance examples: the Cauchy distribution I discussed last time and the Student’s t-distribution with two degrees of freedom that Efrique suggested. Here, the same characterization (mean, median, or mode) is summarized for both distributions in side-by-side boxplots to facilitate comparisons. It is clear from these boxplots that the results for the median and the mode are essentially identical for these distributions, but the results for the mean differ dramatically (recall that these results are truncated for the Cauchy distribution: 13.6% of the 1000 computed means fell outside the +/- 5 range shown here, exhibiting values approaching +/- 1000). This difference illustrates Efrique’s further point that the mean of the data values is a consistent estimator of the (well-defined) population mean of the Student’s t-distribution with 2 degrees of freedom, while it is not a consistent estimator for the Cauchy distribution. Still, it also clear from this plot that the mean is substantially more variable for the Student’s t-distribution with 2 degrees of freedom than either the median or the modeest mode estimate.

Another example of an infinite-variance distribution where the mean is well-defined but highly variable is the Pareto type I distribution, discussed in Section 4.5.8 of Exploring Data. My favorite reference on distributions is the two volume set by Johnson, Kotz, and Balakrishnan (Continuous Univariate Distributions, Vol. 1 (Wiley Series in Probability and Statistics) and Continuous Univariate Distributions, Vol. 2 (Wiley Series in Probability and Statistics)), who devote an entire 55 page chapter (Chapter 20 in Volume 1) to the Pareto distribution, noting that it is named after Vilafredo Pareto, a mid nineteenth- to early twentieth-century Swiss professor of economics, who proposed it as a description of the distribution of income over a population. In fact, there are several different distributions named after Pareto, but the type I distribution considered here exhibits a power-law decay like the Student’s t-distributions, but it is a J-shaped distribution whose mode is equal to its minimum value. More specifically, this distribution is defined by a location parameter that determines this minimum value and a shape parameter that determines how rapidly the tail decays for values larger than this minimum. The example considered here takes this minimum value as 1 and the shape parameter as 1.5, giving a distribution with a finite mean but an infinite variance. As in the above example, the boxplot summary shown below characterizes the mean, median, and mode for 1000 statistically independent random samples drawn from this distribution, each of size N = 100. As before, it is clear from this plot that the mean is much more highly variable than either the median or the mode.

In this case, however, we have the added complication that since this distribution is not symmetric, its mean, median and mode do not coincide. In fact, the population mode is the minimum value (which is 1 here), corresponding to the solid line at the bottom of the plot. The narrow range of the boxplot values around this correct value suggest that the modeest package is reliably estimating this mode value, but as I noted in my last post, this characterization is not useful here because it tells us nothing about the rate at which the density decays. The theoretical median value can also be calculated easily for this distribution, and here it is approximately equal to 1.587, corresponding to the dashed horizontal line in the plot. As with the mode, it is clear from the boxplot that the median estimated from the data is in generally excellent agreement with this value. Finally, the mean value for this particular distribution is 3, corresponding to the dotted horizontal line in the plot. Since this line lies fairly close to the upper quartile of the computed means (i.e., the top of the “box” in the boxplot), it follows that the estimated mean falls below the correct value almost 75% of the time, but it is also clear that when the mean is overestimated, the extent of this overestimation can be very large. Motivated in part by the fact that the mean doesn’t always exist for the Pareto distribution, Johnson, Kotz and Balakrishnan note in their chapter on these distributions that alternative location measures have been considered, including both the geometric and harmonic means. I will examine these ideas further in a future post.

Finally, klr mentioned my post on useless averages in his blog TimelyPortfolio, where he discusses alternatives to the moving average in characterizing financial time-series. For the case he considers, klr compares a 10-month moving average, the corresponding moving median, and a number of the corresponding mode estimators from the modeest package. This is a very interesting avenue of exploration for me since it is closely related to the median filter and other nonlinear digital filters that can be very useful in cleaning noisy time-series data. I discuss a number of these ideas – including moving-window extensions of other data characterizations like skewness and kurtosis – in my book Mining Imperfect Data: Dealing with Contamination and Incomplete Records.

Again, thanks to all of you for your comments. You have given me much to think about and investigate further, which is one of the joys of doing this blog.

When are averages useless?

2011-08-20T08:21:00.000-07:00

Of all possible single-number characterizations of a data sequence, the average is probably the best known. It is also easy to compute and in favorable cases, it provides a useful characterization of “the typical value” of a sequence of numbers. It is not the only such “typical value,” however, nor is it always the most useful one: two other candidates – location estimators in statistical terminology – are the median and the mode, both of which are discussed in detail in Section 4.1.2 of Exploring Data in Engineering, the Sciences, and Medicine. Like the average, these alternative location estimators are not always “fully representative,” but they do represent viable alternatives – at least sometimes – in cases where the average is sufficiently non-representative as to be effectively useless. As the title of this post suggests, the focus here is on those cases where the mean doesn’t really tell us what we want to know about a data sequence, briefly examining why this happens and what we can do about it.

First, it is worth saying a few words about the two alternatives just mentioned: the median and the mode. Of these, the mode is both the more difficult to estimate and the less broadly useful. Essentially, “the mode” corresponds to “the location of the peak in the data distribution.” One difficulty with this somewhat loose definition is that “the mode” is not always well-defined. The above collection of plots shows three examples where the mode is not well-defined, and another where the mode is well-defined but not particularly useful. The upper left plot shows the density of the uniform distribution on the range [1,2]: there, the density is constant over the entire range, so there is no single, well-defined “peak” or unique maximum to serve as a mode for this distribution. The upper right plot shows a nonparametric density estimate for the Old Faithful geyser waiting time data that I have discussed in several of my recent posts (the R data object faithful). Here, the difficulty is that there are not one but two modes, so “the mode” is not well-defined here, either: we must discuss “the modes.” The same behavior is observed for the arcsin distribution, whose density is shown in the lower left plot in the above figure. This density corresponds to the beta distribution with shape parameters both equal to ½, giving a bimodal distribution whose cumulative probability function can be written simply in terms of the arcsin function, motivating its name (see Section 4.5.1 of Exploring Data for a more complete discussion of both the beta distribution family and the special case of the arcsin distribution). In this case, the two modes of the distribution occur at the extremes of the data, at x = 1 and x = 2.

The second difficulty with the mode noted above is that it is sometimes well-defined but not particularly useful. The case of the J-shaped exponential density shown in the lower right plot above illustrates this point: this distribution exhibits a single, well-defined peak at the minimum value x = 0. Here, you don’t even have to look at the data to arrive at this result, which therefore tells you nothing about the data distribution: this density is described by a single parameter that determines how slowly or rapidly the distribution decays and the mode is independent of this parameter. Despite these limitations, there are cases where the mode represents an extremely useful data characterization, even though it is much harder to estimate than the mean or the median. Fortunately, there is a nice package available in R to address this problem: the modeest package provides 11 different mode estimation procedures. I will illustrate one of these in the examples that follow – the half range mode estimator of Bickel – and I will give a more complete discussion of this package in a later post.

The median is a far better-known data characterization than the mode, and it is both much easier to estimate and much more broadly applicable. In particular, unlike either the mean or the mode, the median is well-defined for any proper data distribution, a result demonstrated in Section 4.1.2 of Exploring Data. Conceptually, computing the median only requires sorting the N data values from smallest to largest and then taking either the middle element from this sorted list (if N is odd), or averaging the middle two elements (if N is even).

The mean is, of course, both the easiest of these characterizations to compute – simply add the N data values and divide by N – and unquestionably the best known. There are, however, at least three situations where the mean can be so highly non-representative as to be useless:

1.      if severe outliers are present;
2.      if the distribution is multi-modal;
3.      if the distribution has infinite variance.

The rest of this post examines each of these cases in turn.

I have discussed the problem of outliers before, but they are an important enough problem in practice to bear repeating. (I devote all of Chapter 7 to this topic in Exploring Data.) The plot below shows the makeup flow rate dataset, available from the companion website for Exploring Data (the dataset is makeup.csv, available on the R programs and datasets page). This dataset consists of 2,589 successive measurements of the flow rate of a fluid stream in an industrial manufacturing process. The points in this plot show two distinct forms of behavior: those with values on the order of 400 represent measurements made during normal process operation, while those with values less than about 300 correspond to measurements made when the process is shut down (these values are approximately zero) or is in the process of being either shut down or started back up. The three lines in this plot correspond to the mean (the solid line at approximately 315), the median (the dotted line at approximately 393), and the mode (the dashed line at approximately 403, estimated using the “hrm” method in the modeest package). As I have noted previously, the mean in this case represents a useful line of demarcation between the normal operation data (those points above the mean, representing 77.6% of the data) and the shutdown segments (those points below the mean, representing 22.4% of the data). In contrast, both the median and the specific mode estimator used here provide much better characterizations of the normal operating data.

The next plot below shows a nonparametric density estimate of the Old Faithful geyser waiting data I discussed in my last few posts. The solid vertical line at 70.90 corresponds to the mean value computed from the complete dataset. It has been said that a true compromise is an agreement that makes all parties equally unhappy, and this seems a reasonable description of the mean here: the value lies about mid-way between the two peaks in this distribution, centered at approximately 55 and 80; in fact, this value lies fairly close to the trough between the peaks in this density estimate. (The situation is even worse for the arcsin density discussed above: there, the two modes occur at values of 1 and 2, while the mean falls equidistant from both at 1.5, arguably the “least representative” value in the whole data range.) The median waiting time value is 76, corresponding to the dotted line just to the left of the main peak at about 80, and the mode (again, computed using the package modeest with the “hrm” method) corresponds to the dashed line at 83, just to the right of the main peak. The basic difficulty here is that all of these location estimators are inherently inadequate since they are attempting to characterize “the representative value” of a data sequence that has “two representative values:” one representing the smaller peak at around 55 and the other representing the larger peak at around 80. In this case, both the median and the mode do a better job of characterizing the larger of the two peaks in the distribution (but not a great job), although such a partial characterization is not always what we want. This type of behavior is exactly what the mixture models I discussed in my last few posts are intended to describe.

To illustrate the third situation where the mean is essentially useless, consider the Cauchy distribution, corresponding to the Student’s t distribution with one degree of freedom. This is probably the best known infinite-variance distribution there is, and it is often used as an extreme example because it causes a lot of estimation procedures to fail. The plot below is a (truncated) boxplot comparison of the values of the mean, median, and mode computed from 1000 independently generated Cauchy random number sequences, each of length N = 100. It is clear from these boxplots that the variability of the mean is much greater than that of either of the other two estimators, which are the median and the mode, the latter again estimated from the data using the half-range mode (hrm) method in the modeest package. One of the consequences of working with infinite variance distributions is that the mean is no longer a consistent location estimator, meaning that the variance of the estimated mean does not approach zero in the limit of large sample sizes. In fact, the Cauchy distribution is one of the examples I discuss in Chapter 6 of Exploring Data as a counterexample to the Central Limit Theorem: for most data distributions, the distribution of the mean approaches a Gaussian limit with a variance that decreases inversely with the sample size N, but for the Cauchy distribution, the distribution of the mean is exactly the same as that of the data itself. In other words, for the Cauchy distribution, averaging a collection of N numbers does not reduce the variability at all. This is exactly what we are seeing here, although the plot below doesn’t show how bad the situation really is: the smallest value of the mean in this sequence of 1000 estimates is -798.97 and the largest value is 928.85. In order to see any detail at all in the distribution of the median and mode values, it was necessary to restrict the range of the boxplots shown here to lie between -5 and +5, which eliminated 13.6% of the computed mean values. In contrast, the median is known to be a reasonably good location estimator for the Cauchy distribution (see Section 6.6.1 of Exploring Data for a further discussion of this point), and the results presented here suggest that Bickel’s half-range mode estimator is also a reasonable candidate. The main point here is that the mean is a completely unreasonable estimator in situations like this one, an important point in view of the growing interest in data models like the infinite-variance Zipf distribution to describe “long-tailed” phenomena in business.

I will have more to say about both the modeest package and Zipf distributions in upcoming posts.

Fitting mixture distributions with the R package mixtools

2011-08-06T14:23:00.000-07:00

My last two posts have been about mixture models, with examples to illustrate what they are and how they can be useful. Further discussion and more examples can be found in Chapter 10 of Exploring Data in Engineering, the Sciences, and Medicine. One important topic I haven’t covered is how to fit mixture models to datasets like the Old Faithful geyser data that I have discussed previously: a nonparametric density plot gives fairly compelling evidence for a bimodal distribution, but how do you estimate the parameters of a mixture model that describes these two modes? For a finite Gaussian mixture distribution, one way is by trial and error, first estimating the centers of the peaks by eye in the density plot (these become the component means), and adjusting the standard deviations and mixing percentages to approximately match the peak widths and heights, respectively. This post considers the more systematic alternative of estimating the mixture distribution parameters using the mixtools package in R.

The mixtools package is one of several available in R to fit mixture distributions or to solve the closely related problem of model-based clustering. Further, mixtools includes a variety of procedures for fitting mixture models of different types. This post focuses on one of these – the normalmixEM procedure for fitting normal mixture densities – and applies it to two simple examples, starting with the Old Faithful dataset mentioned above. A much more complete and thorough discussion of the mixtools package – which also discusses its application to the Old Faithful dataset – is given in the R package vignette, mixtools: An R Package for Analyzing Finite Mixture Models.

The above plot shows the results obtained using the normalmixEM procedure with its default parameter values, applied to the Old Faithful waiting time data. Specifically, this plot was generated by the following sequence of R commands:

            library(mixtools)
            wait = faithful$waiting
            mixmdl = normalmixEM(wait)
            plot(mixmdl,which=2)
            lines(density(wait), lty=2, lwd=2)

Like many modeling tools in R, the normalmixEM procedure has associated plot and summary methods. In this case, the plot method displays either the log likelihood associated with each iteration of the EM fitting algorithm (more about that below), or the component densities shown above, or both. Specifying “which=1” displays only the log likelihood plot (this is the default), specifying “which = 2” displays only the density components/histogram plot shown here, and specifying “density = TRUE” without specifying the “which” parameter gives both plots. Note that the two solid curves shown in the above plot correspond to the individual Gaussian density components in the mixture distribution, each scaled by the estimated probability of an observation being drawn from that component distribution. The final line of R code above overlays the nonparametric density estimate generated by the density function with its default parameters, shown here as the heavy dashed line (obtained by specifying “lty = 2”).

Most of the procedures in the mixtools package are based on the iterative expectation maximization (EM) algorithm, discussed in Section 2 of the mixtools vignette and also in Chapter 16 of Exploring Data. A detailed discussion of this algorithm is beyond the scope of this post – books have been devoted to the topic (see, for example, the book by McLachlan and Krishnan, The EM Algorithm and Extensions (Wiley Series in Probability and Statistics) ) – but the following two points are important to note here. First, the EM algorithm is an iterative procedure, and the time required for it to reach convergence – if it converges at all – depends strongly on the problem to which it is applied. The second key point is that because it is an iterative procedure, the EM algorithm requires starting values for the parameters, and algorithm performance can depend strongly on these initial values. The normalmixEM procedure supports both user-supplied starting values and built-in estimation of starting values if none are supplied. These built-in estimates are the default and, in favorable cases, they work quite well. The Old Faithful waiting time data is a case in point – using the default starting values gives the following parameter estimates:

            > mixmdl[c("lambda","mu","sigma")]
$lambda
[1] 0.3608868 0.6391132

$mu
[1] 54.61489 80.09109

$sigma
[1] 5.871241 5.867718

The mixture density described by these parameters is given by:

p(x) = lambda[1] n(x; mu[1], sigma[1]) + lambda[2] n(x; mu[2], sigma[2])

where n(x; mu, sigma) represents the Gaussian probability density function with mean mu and standard deviation sigma.

One reason the default starting values work well for the Old Faithful waiting time data is that if nothing is specified, the number of components (the parameter k) is set equal to 2. Thus, if you are attempting to fit a mixture model with more than two components, this number should be specified, either by setting k to some other value and not specifying any starting estimates for the parameters lambda, mu, and sigma, or by specifying a vector with k components as starting values for at least one of these parameters. (There are a number of useful options in calling the normalmixEM procedure: for example, specifying the initial sigma value as a scalar constant rather than a vector with k components forces the component variances to be equal. I won’t attempt to give a detailed discussion of these options here; for that, type “help(normalmixEM)”.)

Another important point about the default starting values is that, aside from the number of components k, any unspecified initial parameter estimates are selected randomly by the normalmixEM procedure. This means that, even in cases where the default starting values consistently work well – again, the Old Faithful waiting time dataset seems to be such a case – the number of iterations required to obtain the final result can vary significantly from one run to the next. (Specifically, the normalmixEM procedure does not fix the seed for the random number generators used to compute these starting values, so repeated runs of the procedure with the same data will start from different initial parameter values and require different numbers of iterations to achieve convergence. In the case of the Old Faithful waiting time data, I have seen anywhere between 16 and 59 iterations required, with the final results differing only very slightly, typically in the fifth or sixth decimal place. If you want to use the same starting value on successive runs, this can be done by setting the random number seed via the set.seed command before you invoke the normalmixEM procedure.)

It is important to note that the default starting values do not always work well, even if the correct number of components is specified. This point is illustrated nicely by the following example. The plot above shows two curves: the solid line is the exact density for the three-component Gaussian mixture distribution described by the following parameters:

            mu = (2.00, 5.00, 7.00)
            sigma = (1.000, 1.000, 1.000)
            lambda = (0.200, 0.600, 0.200)

The dashed curve in the figure is the nonparametric density estimate generated from n = 500 observations drawn from this mixture distribution. Note that the first two components of this mixture distribution are evident in both of these plots, from the density peaks at approximately 2 and 5. The third component, however, is too close to the second to yield a clear peak in either density, giving rise instead to slightly asymmetric “shoulders” on the right side of the upper peaks. The key point is that the components in this mixture distribution are difficult to distinguish from either of these density estimates, and this hints at further difficulties to come.

Applying the normalmixEM procedure to the 500 sample sequence used to generate the nonparametric density estimate shown above and specifying k = 3 gives results that are substantially more variable than the Old Faithful results discussed above. In fact, to compare these results, it is necessary to be explicit about the values of the random seeds used to initialize the parameter estimation procedure. Specifying this random seed as 101 and only specifying k=3 in the normalmixEM call yields the following parameter estimates after 78 iterations:

            mu = (1.77, 4.87, 5.44)
            sigma = (0.766, 0.115, 1.463)
            lambda = (0.168, 0.028, 0.803)

Comparing these results with the correct parameter values listed above, it is clear that some of these estimation errors are quite large. The figure shown below compares the mixture density constructed from these parameters (the heavy dashed curve) with the nonparametric density estimate computed from the data used to estimate them. The prominent “spike” in this mixture density plot corresponds to the very small standard deviation estimated for the second component and it provides a dramatic illustration of the relatively poor results obtained for this particular example.

Repeating this numerical experiment with different random seeds to obtain different random starting estimates, the normalmixEM procedure failed to converge in 1000 iterations for seed values of 102 and 103, but it converged after 393 iterations for the seed value 104, yielding the following parameter estimates:

            mu = (1.79, 5.03, 5.46)
            sigma = (0.775, 0.352, 1.493)
            lambda = (0.169, 0.063, 0.768)

Arguably, the general behavior of these parameter estimates is quite similar to those obtained with the random seed value 101, but note that the second variance component differs by a factor of three, and the second component of lambda increases almost as much.

Increasing the sample size from n = 500 to n = 2000 and repeating these experiments, the normalmixEM procedure failed to converge after 1000 iterations for all four of the random seed values 101 through 104. If, however, we specify the correct standard deviations (i.e., specify “sigma = c(1,1,1)” when we invoke normalmixEM) and we increase the maximum number of iterations to 3000 (i.e., specify “maxit = 3000”), the procedure does converge after 2417 iterations for the seed value 101, yielding the following parameter estimates:

            mu = (1.98, 4.98, 7.15)
            sigma = (1.012, 1.055, 0.929)
            lambda = (0.198, 0.641, 0.161)

While these parameters took a lot more effort to obtain, they are clearly much closer to the correct values, emphasizing the point that when we are fitting a model to data, our results generally improve as the amount of available data increases and as our starting estimates become more accurate. This point is further illustrated by the plot shown below, analogous to the previous one, but constructed from the model fit to the longer data sequence and incorporating better initial parameter estimates. Interestingly, re-running the same procedure but taking the correct means as starting parameter estimates instead of the correct standard deviations, the procedure failed to converge in 3000 iterations.

Overall, I like what I have seen so far of the mixtools package, and I look forward to exploring its capabilities further. It’s great to have a built-in procedure – i.e., one I didn’t have to write and debug myself – that does all of the things that this package does. However, the three-component mixture results presented here do illustrate an important point: the behavior of iterative procedures like normalmixEM and others in the mixtools package can depend strongly on the starting values chosen to initialize the iteration process, and the extent of this dependence can vary greatly from one application to another.

Mixture distributions and models: a clarification

2011-07-16T11:32:00.000-07:00

In response to my last post, Chris had the following comment:

I am actually trying to better understand the distinction between mixture models and mixture distributions in my own work. You seem to say mixture models apply to a small set of models – namely regression models.

This comment suggests that my caution about the difference between mixed-effect models and mixture distributions may have caused as much confusion as clarification, and the purpose of this post is to try to clear up this confusion.

So first, let me offer the following general observations. The terms “mixture models” refers to a generalization of the class of finite mixture distributions that I discussed in my previous post. I give a more detailed discussion of finite mixture distributions in Chapter 10 of Exploring Data in Engineering, the Sciences, and Medicine , and the more general class of mixture models is discussed in the book Mixture Models (Statistics: A Series of Textbooks and Monographs) by Geoffrey J. McLachlan and Kaye E. Bashford. The basic idea is that we are describing some observed phenomenon like the Old Faithful geyser data (the faithful data object in R) where a close look at the data (e.g., with a nonparametric density estimate) suggests substantial heterogeneity. In particular, the density estimates I presented last time for both of the variables in this dataset exhibit clear evidence of bimodality. Essentially, the idea behind a mixture model/mixture distribution is that we are observing something that isn’t fully characterized by a single, simple distribution or model, but instead by several such distributions or models, with some random selection mechanism at work. In the case of mixture distributions, some observations appear to be drawn from distribution 1, some from distribution 2, and so forth. The more general class of mixture models is quite broad, including things like heterogeneous regression models, where the response may depend approximately linearly on some covariate with one slope and intercept for observations drawn from one sub-population, but with another, very different slope and intercept for observations drawn from another sub-population. I present an example at the end of this post that illustrates this idea.

The probable source of confusion for Chris – and very possibly other readers – is the comment I made about the difference between these mixture models and mixed-effect models. This other class of models – which I only mentioned in passing in my post – typically consists of a linear regression model with two types of prediction variables: deterministic predictors, like those that appear in standard linear regression models, and random predictors that are typically assumed to obey a Gaussian distribution. This framework has been extended to more general settings like generalized linear models (e.g., mixed-effect logistic regression models). The R package lme4 provides support for fitting both linear mixed-effect models and generalized linear mixed-effect models to data. As I noted last time, these model classes are distinct from the mixture distribution/mixture model classes I discuss here. The models that I do discuss – mixture models – have strong connections with cluster analysis, where we are given a heterogeneous group of objects and typically wish to determine how many distinct groups of objects are present and assign individuals to the appropriate groups. A very high-level view of the many R packages available for clustering – some based on mixture model ideas and some not – is available from the CRAN clustering task view page. Two packages from this task view that I plan to discuss in future posts are flexmix and mixtools, both of which support a variety of mixture model applications. The following comments from the vignette FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R give an indication of the range of areas where these ideas are useful:

“Finite mixture models have been used for more than 100 years, but have seen a real boost in popularity over the last decade due to the tremendous increase in available computing power. The areas of application of mixture models range from biology and medicine to physics, economics, and marketing. On the one hand, these models can be applied to data where observations originate from various groups and the group affiliations are not known, and on the other hand to provide approximations for multi-modal distributions.”

The following example illustrates the second of these ideas, motivated by the Old Faithful geyser data that I discussed last time. As a reminder, the plot above shows the nonparametric density estimate generated from the 272 observations of the Old Faithful waiting time data included in the faithful data object, using the density procedure in R with the default parameter settings. As I noted last time, the plot shows two clear peaks, the lower one centered at approximately 55 minutes, and the second at approximately 80 minutes. Also, note that the first peak is substantially smaller in amplitude and appears to be somewhat narrower than the second peak.

To illustrate the connection with finite mixture distributions, the R procedure described below generates a two-component Gaussian mixture density whose random samples exhibit approximately the same behavior seen in the Old Faithful waiting time data. The results generated by this procedure are shown in the above figure, which includes two overlaid plots: one corresponding to the exact density for the two-component Gaussian mixture distribution (the solid line), and the other corresponding to the nonparametric density estimate computed from N = 272 random samples drawn from this mixture distribution (the dashed line). As in the previous plot, the nonparametric density estimate was computed using the density command in R with its default parameter values. The first component in this mixture has mean 54.5 and standard deviation 8.0, values chosen by trial and error to approximately match the lower peak in the Old Faithful waiting time distribution. The second component has mean 80.0 and standard deviation 5.0, chosen to approximately match the second peak in the waiting time distribution. The probabilities associated with the first and second components are 0.45 and 0.55, respectively, selected to give approximately the same peak heights seen in the waiting time density estimate. Combining these results, the density of this mixture distribution is:

p(x) = 0.45 n(x; 54.5, 8.0) + 0.55 n(x; 80.0, 5.0),

where n(x;m,s) denotes the Gaussian density function with mean m and standard deviation s. These density functions can be generated using the dnorm function in R.

The R procedure listed below generates n independent, identically distributed random samples from an m-component Gaussian mixture distribution. This procedure is called with the following parameters:

n = the number of random samples to generate

mvec = vector of m mean values

svec = vector of m standard deviations

pvec = vector of probabilities for each of the m components

iseed = integer seed to initialize the random number generators

The R code for the procedure looks like this:

MixEx01GenProc <- function(n, muvec, sigvec, pvec, iseed=101){

set.seed(iseed)

m <- length(pvec)

indx <- sample(seq(1,m,1), size=n, replace=T, prob=pvec)

yvec <- 0

for (i in 1:m){

xvec <- rnorm(n, mean=muvec[i], sd=sigvec[i])

yvec <- yvec + xvec * as.numeric(indx == i)

}

yvec

}

The first statement initializes the random number generator using the iseed parameter, which is given a default value of 101. The second line determines the number of components in the mixture density from the length of the pvec parameter vector, and the third line generates a random sequence indx of component indices taking the values 1 through m with probabilities determined by the pvec parameter. The rest of the program is a short loop that generates each component in turn, using indx to randomly select observations from each of these components with the appropriate probability. To see how this works, note that the first pass through the loop generates the random vector xvec of length n, with mean given by the first element of the vector muvec and standard deviation given by the first element of the vector sigvec. Then, for every one of the n elements of yvec for which the indx vector is equal to 1, yvec is set equal to the corresponding element of this first random component xvec. On the second pass through the loop, the second random component is generated as xvec, again with length n but now with mean specified by the second element of muvec and standard deviation determined by the second element of sigvec. As before, this value is added to the initial value of yvec whenever the selection index vector indx is equal to 2. Note that since every element of the indx vector is unique, none of the nonzero elements of yvec computed during the first iteration of the loop are modified; instead, the only elements of yvec that are modified in the second pass through the loop have their initial value of zero, specified in the line above the start of the loop. More generally, each pass through the loop generates the next component of the mixture distribution and fills in the corresponding elements of yvec as determined by the random selection index vector indx.

As I noted at the beginning of this post, the notion of a mixture model is more general than that of the finite mixture distributions just described, but closely related. I conclude this post with a simple example of a more general mixture model. The above scatter plot shows two variables, x and y, related by the following mixture model:

y = x + e₁ with probability p₁ = 0.40,

and

y = -x + 2 + e₂ with probability p₂ = 0.60,

where e₁ is a zero-mean Gaussian random variable with standard deviation 0.1, and e₂ is a zero-mean Gaussian random variable with standard deviation 0.3. To emphasize the components in the mixture model, points corresponding to the first component are plotted as solid circles, while points corresponding to the second component are plotted as open triangles. The two dashed lines in this plot represent the ordnary least squares regression lines fit to each component separately, and they both correspond reasonably well to the underlying linear relationships that define the two components (e.g., the least squares line fit to the solid circles has a slope of approximately +1 and an intercept of approximately 0). In contrast, the heavier dotted line represents the ordinary least squares regression line fit to the complete dataset without any knowledge of its underlying component structure: this line is almost horizontal and represents a very poor approximation to the behavior of the dataset.

The point of this example is to illustrate two things. First, it provides a relatively simple illustration of how the mixture density idea discussed above generalizes to the setting of regression models and beyond: we can construct fairly general mixture models by requiring different randomly selected subsets of the data to conform to different modeling assumptions. The second point – emphasized by the strong disagreement between the overall regression line and both of the component regression lines – is that if we are given only the dataset (i.e., the x and y values themselves) without knowing which component they represent, standard analysis procedures are likely to perform very badly. This question – how do we analyze a dataset like this one without detailed prior knowledge of its heterogeneous structure – is what R packages like flexmix and mixtools are designed to address.

A Brief Introduction to Mixture Distributions

2011-06-18T11:20:00.000-07:00

Last time, I discussed some of the advantages and disadvantages of robust estimators like the median and the MADM scale estimator, noting that certain types of datasets – like the rainfall dataset discussed last time – can cause these estimators to fail spectacularly. An extremely useful idea in working with datasets like this one is that of mixture distributions, which describe random variables that are drawn from more than one parent population. I discuss mixture distributions in some detail in Chapter 10 of Exploring Data in Engineering, the Sciences, and Medicine (Section 10.6), and the objective of this post is to give a brief introduction to the basic ideas of mixture distributions, with some hints about how they can be useful in connection with problems like analyzing the rainfall data discussed last time. Subsequent posts will discuss these ideas in more detail, with pointers to specific R packages that are useful in applying these ideas to real data analysis problems.

Before proceeding, it is important to emphasize that the topic of this post and subsequent discussions is “mixture distributions” and not “mixture models,” which are the subject of considerable current interest in the statistics community. The distinction is important because these topics are very different: mixture distributions represent a useful way of describing heterogeneity in the distribution of a variable, whereas mixture models provide a foundation for incorporating both deterministic and random predictor variables in regression models.

To motivate the ideas presented here, the figure above shows four plots constructed from the Old Faithful geyser data that I have discussed previously (the R dataset faithful). This data frame contains 272 observations of the durations of successive eruptions and the waiting time until the next eruption, both measured in minutes. The upper left plot above is the normal Q-Q plot for the duration values, computed using the qqPlot procedure from the car package that I have discussed previously, and the upper right plot is the corresponding normal Q-Q plot for the waiting times. Both of these Q-Q plots exhibit pronounced “kinks,” which I have noted previously often indicate the presence of a multimodal data distribution. The lower two plots are the corresponding nonparametric density estimates (computed using the density procedure in R with its default parameters). The point of these plots is that they give strong evidence that the distributions of both the durations and waiting times are bimodal.

In fact, bimodal and more general multimodal distributions arise frequently in practice, particularly in cases where we are observing a composite response from multiple, distinct sources. This is the basic idea behind mixture distributions: the response x that we observe is modeled as a random variable that has some probability p₁ of being drawn from distribution D₁, probability p₂ of being drawn from distribution D₂, and so forth, with probability p_n of being drawn from distribution D_n, where n is the number of components in our mixture distribution. The key assumption here is one of statistical independence between the process of randomly selecting the component distribution D_i to be drawn and these distributions themselves. That is, we assume there is a random selection process that first generates the numbers 1 through n with probabilities p₁ through p_n. Then, once we have drawn some number j, we turn to distribution D_j and draw the random variable x from this distribution. So long as the probabilities – or mixing percentages – p_i sum to 1, and all of the distributions D_i are proper densities, the combination also defines a proper probability density function, which can be used as the basis for computing expectations, formulating maximum likelihood estimation problems, and so forth. The simplest case is that for n = 2, which provides many different specific examples that have been found to be extremely useful in practice. Because it is the easiest to understand, this post will focus on this special case.

Since the mixing percentages must sum to 1, it follows that p₂ = 1 – p₁ when n = 2, it simplifies the discussion to drop the subscripts, writing p₁ = p and p₂ = 1 – p. Similarly, let f(x) denote the density associated with the distribution D₁ and g(x) denote the density associated with the distribution D₂. The overall probability density function p(x) is then given by:

p(x) = p f(x) + (1 – p)g(x)

To illustrate the flexibility of this idea, consider the case where both f(x) and g(x) are Gaussian distributions. The figure below shows four specific examples.

The upper left plot shows the standard normal distribution (i.e., mean 0 and standard deviation 1), which corresponds to taking both f(x) and g(x) as standard normal densities, for any choice of p. I have included it here both because it represents what has been historically the most common distributional assumption, and because it represents a reference case in interpreting the other distributions considered here. The upper right plot corresponds to the contaminated normal outlier distribution, widely adopted in the robust statistics literature. There, the idea is that measurement errors – traditionally modeled as zero-mean Gaussian random variables with some unknown standard deviation S – mostly conform to this model, but some fraction (typically, 10% to 20%) of these measurement errors have larger variability. The traditional contaminated normal model defines p as this contamination percentage and assumes a standard deviation of 3S for these measurements. The upper right plot above shows the density for this contaminated normal model with p = 0.15 (i.e., 15% contamination). Visually, this plot looks identical to the one on the upper left for the standard normal distribution: plotting them on common axes (as done in Fig. 10.17 of Exploring Data) shows that these distributions are not identical, a point discussed further below on the basis of Q-Q plots.

The lower left plot in the figure above corresponds to a two-component Gaussian mixture distribution where both components are equally represented (i.e., p = 0.5), the first component f(x) is the standard normal distribution as before, and the second component g(x) is a Gaussian distribution with mean 3 and standard deviation 3. The first component contributes the sharper main peak centered at zero, while the second component contributes the broad “shoulder” seen in the right half of this plot. Finally, the lower right plot shows a mixture distribution with p = 0.40 where the first component is a Gaussian distribution with mean -2 and standard deviation 1, and the second component is a Gaussian distribution with mean +2 and standard deviation 1. The result is a bimodal distribution with the same general characteristics as the Old Faithful geyser data.

The point of these examples has been to illustrate the flexibility of the mixture distribution concept, in describing everything from outliers to the natural heterogeneity of natural phenomena with more than one distinct generation mechanism. Before leaving this discussion, it is worth considering the contaminated normal case a bit further. The plot above shows the normal Q-Q plot for 272 samples of the contaminated normal distribution just described, generated using the R procedure listed below (the number 272 was chosen to make the sample size the same as that of the Old Faithful geyser data discussed above). As before, this Q-Q plot was generated using the procedure qqPlot from the car package. The heavy-tailed, non-Gaussian character of this distribution is evident from the fact that both the upper and lower points in this plot fall well outside the 95% confidence interval around the Gaussian reference line. This example illustrates the power of the Q-Q plot for distributional assessment: like the density plots shown above for the standard normal and contaminated normal distributions, nonparametric density plots (not shown) generated from 272 samples drawn from each of these data distributions are not markedly different.

The contaminated normal random samples used to construct this Q-Q plot were generated with the following simple R function:

cngen.proc01 <- function(n=272,cpct = 0.15, mu1 = 0, mu2 = 0, sig1 = 1, sig2 = 3,iseed=101){

set.seed(iseed)

y0 <- rnorm(n,mean=mu1, sd = sig1)

y1 <- rnorm(n,mean=mu2, sd = sig2)

flag <- rbinom(n,size=1,prob=cpct)

y <- y0*(1 - flag) + y1*flag

}

This function is called with seven parameters, all of which are given default values. The parameter n is the sample size, cpct is the contamination percentage (expressed as a fraction, so cpct = 0.15 corresponds to 15% contamination), mu1 and mu2 are the means of the component Gaussian distributions, and sig1 and sig2 are the corresponding standard deviations. The default values for cpct, mu1, mu2, sig1, and sig2 are those for a typical contaminated normal outlier distribution. Finally, the parameter iseed is the seed for the random number generator: specifying its value as a default means that the procedure will return the same set of pseudorandom numbers each time it is called; to obtain an independent set, simply specify a different value for iseed.

The procedure itself generates three mutually independent random variables: y0 corresponds to the first component distribution, y1 corresponds to the second, and flag is a binomial random variable that determines which component distribution is selected for the final random number: when flag = 1 (an event with probability cpct), the contaminated value y1 is returned, and when flag = 0 (an event with probability 1 – cpct), the non-contaminated value y0 is returned.

The basic ideas just described – random selection of a random variable from n distinct distributions – can be applied to a very wide range of distributions, leading to an extremely wide range of data models. The examples described above give a very preliminary indication of the power of Gaussian mixture distributions as approximations to the distribution of heterogeneous phenomena that arise from multiple sources. Practical applications of more general (i.e., not necessarily Gaussian) mixture distributions include modeling particle sizes generated by multiple mechanisms (e.g., accretion of large particles from smaller ones and fragmentation of larger particles to form smaller ones, possibly due to differences in material characteristics), pore size distributions in rocks, polymer chain length or chain branching distributions, and many others. More generally, these ideas are also applicable to discrete distributions, leading to ideas like the zero-inflated Poisson and negative binomial distributions that I will discuss in my next post, or to combinations of continuous distributions with degenerate distributions (e.g., concentrated at zero), leading to ideas like zero-augmented continuous distributions that may be appropriate in applications like the analysis of the rainfall data I discussed last time. One of the advantages of R is that it provides a wide variety of support for various useful implementations of these ideas.

The pros and cons of robust data characterizations

2011-06-06T17:40:00.000-07:00

Over the years, I have looked at a lot of data contaminated with outliers, the subject of Chapter 7 of Exploring Data in Engineering, the Sciences, and Medicine. That chapter adopts the definition of an outlier presented by Barnett and Lewis in their book Outliers in Statistical Data 2nd Edition, that outliers are “data points inconsistent with the majority of values in a dataset.” The practical importance of outliers lies in the fact that even a few of them can cause standard data characterization and analysis procedures to give highly misleading results. Robust procedures have been developed to address these problems, and their performance in the face of outliers can be dramatically better than that of standard methods. On the other hand, robust procedures can also fail, and when they do, they often fail spectacularly.

The above figure is a plot of successive measurements of the makeup flow rate from an industrial manufacturing process. This example is introduced in Chapter 2 of Exploring Data and is discussed in detail in Chapter 7; the dataset itself is available as the dataset makeup.csv from the companion website

for the book. Each point in this dataset corresponds to a measurement of the flow rate of a liquid component in a process with a recycle stream: most of this component is recovered from a downstream processing step and recycled back to the upstream step, but a small amount of this material is lost to evaporation and must be “made up” by a separate feed stream. The data points plotted above correspond to the flow rate of this makeup stream. The quantile function in R computes the Tukey five number summary of the data sequence (i.e., the minimum value, lower quartile, median, upper quartile, and maximum value). Applied to this sequence, the results are:

> quantile(MakeUpFlow)

0% 25% 50% 75% 100%

0.00086 370.43747 393.35861 404.29602 439.59091

Looking at the plot, it is clear that these observations partition naturally into two subsets: those values between 0 and about 200 – all falling between the minimum value and the lower quartile – and those values between about 300 and 440. In fact, 12.2% of the data values lie between 0 and 10, and 22.3% of the values are less than 200. Physically, these smaller-than-average flow rates correspond to shutdown episodes, where the manufacturing process is either not running, is in the process of being shut down, or is in the process of being started back up. If we compute the mean of the complete dataset, we obtain an average makeup flow rate of 315.46, corresponding to the dotted line in the above plot. It is clear that while this number represents the center of the overall data distribution, it is not representative of either the shutdown episodes or normal process operation, since it forms a perfect dividing line between these two data subsets. If we exclude the shutdown episodes – specifically, if we omit all values smaller than 200 – and compute the mean of what remains, we obtain 397.56, corresponding to the solid line in the plot above. It is clear that this number much better represents “typical process operation.”

This example illustrates the well-known outlier sensitivity of the mean. An outlier-resistant alternative is the median, defined as follows: if the number N of data values is odd, the median is the middle element in the ordered list we obtain by sorting these values from smallest to largest; if N is even, the median is the average of the two middle elements in this list. Applied to the complete makeup flow rate dataset, the median gives a “typical” flow rate of 393.36, fairly close to the average value computed without the shutdown episodes. Excluding these episodes changes the median value only a little, to 399.50, again quite close to the mean value after we remove these anomalous data points.

The sensitivity of the standard deviation to outliers is even worse than that of the mean. For the complete dataset, the computed standard deviation is 155.05, an order of magnitude larger than the value (15.22) computed from the dataset without the shutdown episodes. While it is not as well known as the median, an outlier-resistant alternative to the standard deviation is the MADM scale estimator, constructed as follows. First, compute the median as a reference value and compute the differences between every data point and this reference value. Then, take the absolute values of these differences, giving a set of non-negative numbers that tell us how far each point lies from the median reference value. Next, we compute the median of this sequence of absolute differences, which tells us how far a typical data point lies from the reference value. Because we have used the outlier-resistant median both to compute the reference value and to compute this “typical distance” number, the result is highly outlier-resistant. The only difficulty is that, for approximately normally distributed data with no outliers, this number is consistently smaller than the standard deviation. To obtain an alternative estimator of the standard deviation – i.e., a number that is approximately equal to the standard deviation when we compute it from a large, outlier-free collection of approximately Gaussian data – we need to multiply the number just defined by 1.4826. (A derivation of this result is given in Chapter 7 of Exploring Data.) This bias-corrected scale estimate is available as the built-in function mad in R, and applying it to the makeup flow rate dataset gives a scale estimate of 20.22, about the same magnitude as the standard deviation of the data sequence with the shutdown episodes removed. The MADM scale estimate is not completely immune to outliers, so when it is applied to the dataset without the shutdown episodes, it gives a somewhat smaller scale estimate of 13.40. Still, both of these estimates are much more consistent with the standard deviation of the normal operation data than the uncorrected standard deviation computed from the complete dataset.

The key point of this discussion has been to provide a simple, real-data illustration of both the considerable outlier sensitivity of the “standard” location and scale estimators (i.e., the mean and standard deviation) and the existence and performance of robust alternatives to both (i.e., the median and MADM scale estimator). As I noted at the beginning of this discussion, however, while robust estimators can work extremely well – as in this example – they can also fail spectacularly. In general, robust estimators are designed to protect us against anomalous observations that occur out in the tails of the distribution. In the case of “light-tailed” distributions like the uniform distribution, robust measures like the MADM scale estimator can exhibit surprising behavior. In particular, in the absence of outliers, we expect the standard deviation and the MADM scale estimate to give us roughly comparable results. If the standard deviation is much larger than the MADM scale estimate – as in the makeup flow rate example just discussed – this can usually be taken as an indication of either a naturally heavy-tailed data distribution (e.g., the Student’s t distribution with few degrees of freedom), or, as in this example, the presence of contaminating outliers. In the face of light-tailed data distributions, however, the MADM scale estimate can be substantially larger than the standard deviation, and we need to be aware of this fact if we use the MADM scale estimator as a measure of data spread.

More serious problems arise in the case of coarsely quantized data. In particular, it was once popular in engineering applications to record variables like temperature only to single digit accuracy (e.g., temperature to the nearest tenth of a degree). This practice was motivated both by limited memory available in early process monitoring computers, and by a recognition that this was roughly the limit of measurement accuracy, anyway. This kind of data truncation does degrade traditional characterizations like the standard deviation somewhat, but the consequences for the MADM scale estimate can be catastrophic. The plot above shows a sequence of 100 Gaussian data samples, with mean 0 and standard deviation 0.15, represented two ways. The line represents the original Gaussian data sample, while the solid circles correspond to these data values truncated to a single digit (e.g., the first value in the original data sequence is -0.11658569, which gets truncated to -0.1). Both the standard deviation and the MADM scale estimate are reasonably accurate when computed from the unmodified data sequence: the standard deviation is 0.137, while the MADM scale estimate is 0.133, both within about 10% of the correct value of 0.150. If we compute the standard deviation from the truncated data samples, the value declines to 0.103, reflecting the fact that by truncating the data from the 9 digits generated by the random number generator to 1, we are consistently reducing the magnitude of each number, so now our estimation error is more like 30%. Here, however, the MADM scale estimate fails completely, returning the value 0. The reason for this is that, on truncation, 55% of these data values are zero (note that any value between -0.09 and 0.09 will be truncated to zero, and these data values occur quite commonly in a zero-mean, normally distributed data sequence with standard deviation 0.15). Whenever more than 50% of the data values are equal, we don’t need to bother computing either the median or the MADM scale estimator. In particular, the median of this data sequence is simply equal to this majority value. Similarly, since more than 50% of the observed data values lie a distance of zero from the median, the median of the absolute deviation sequence on which the MADM scale estimate is based is necessarily zero. Thus, the MADM scale for any such sequence is zero, as in this example; this behavior is commonly called implosion.

This coarse quantization example indirectly illustrates one of the key distinctions between continuous and discrete variables: for a continuously distributed random variable, no two observations can have exactly the same value, but for discrete random variables, such ties are characteristic. The coarse quantization example represents something in between, being a discrete approximation of a continuous variable. The distinction is important because it is precisely these ties that were responsible for the complete failure of the MADM scale estimator in the example just described. There are other cases where the underlying data distribution has both continuous and discrete components, and these cases can also cause rank- and order-based data characterizations like the MADM scale estimator to behave very badly. The plot above shows rainfall measurements made every half hour for two months (a total of 732 data points), included in the dataset associated with Chapter 14 (“Half-hourly Precipitation and Streamflow, River Hirnant, Wales, U.K., November and December, 1972”) from the book Data:: A Collection of Problems from Many Fields for the Student and Research Worker (Springer Series in Statistics), by D.F. Andrews and A.M. Herzberg. It was noted above that the built-in function quantile in R generates the Tukey five-number summary for a data sequence, but by specifying the prob parameter, it is possible to obtain finer-grained data characterizations. For example, the following command gives the deciles of this rainfall dataset:

>quantile(Rain,prob=seq(0,1,0.1))

0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.040 0.146 0.440 2.800

It is clear from these results that more than 60% of these recorded rainfall values are identically zero, reflecting the fact that it does not rain most days. As with the coarse quantization example discussed above, the MADM scale estimator implodes for this example, returning the value zero. Similarly, the median value of these numbers is also zero, but the mean is not (specifically, it is 0.136), and the standard deviation is not (this value is 0.354).

Even though the mean and standard deviation give more reasonable-sounding results for the rainfall dataset than their robust counterparts do, it is not completely clear how to interpret these numbers. In particular, we might be interested in asking any or all of the following three types of question:

1. What is the typical rainfall for this region in a month? (or in a week, in the entire two-month period, in a day, etc.)

2. How frequently does it rain?

3. When it does rain, what is the average amount?

Characterizations like the mean and standard deviation may be used to address the first of these three questions, but not the other two. As a preliminary answer to the second question, we can note that rainfall occurs in 34.6% of the half-hourly intervals recorded here, but this requires a slightly more detailed look at the data than what the mean and standard deviation give us. Additional work would be required if we wanted an estimate of how probable it was to rain on any given day, since we would need to summarize the data at the level of days instead of half-hourly observations. Finally, to address questions of the third type, it is necessary to segment the dataset into a “rainy part” and a “dry part.” Note that this segmentation is analogous to that made for the makeup flow rate dataset, but the difference is that here we are interested in characterizing the minority part of the dataset (i.e., “the outliers”) instead of the majority part, which is what the robust estimators are good at. Special types of statistical models have been developed to deal with these situations, and I will discuss some of these in my next post. For now, I will conclude with two key observations.

The first is that of Collin Mallows (C.L. Mallows, “Robust methods – some examples of their use,” American Statistician, vol. 33, 1979, pages 179-184), which I quote at the end of Chapter 7 of Exploring Data:

“A simple and useful strategy is to perform one’s analysis both robustly and by standard methods and to compare the results. If the differences are minor, either set can be presented. If the differences are not, one must perforce consider why not, and the robust analysis is already at hand to guide the next steps.

The importance of these considerations is enhanced when we are dealing with large amounts of data, since then examining all the data in detail is impractical, and we are forced to contemplate working with data that is, almost certainly, partially bad and with models that are almost certainly inadequate.”

The second observation is that Mallows advice is practical for almost any type of data analysis we want to consider, since robust methods exist to handle a very wide variety of different analysis tasks. An excellent introduction to many of these robust methods is given by Rand Wilcox in his book Introduction to Robust Estimation and Hypothesis Testing, Second Edition (Statistical Modeling and Decision Science), which includes R procedures for implementing most of the methods he discusses.

The distribution of interestingness

2011-05-21T10:51:00.000-07:00

On April 22, David Landy posed a question about the distribution of interestingness values in response to my April 3^rd post on “Interestingness Measures.” He noted that the survey paper by Hilderman and Hamilton that I cited there makes the following comment:

“Our belief is that a useful measure of interestingness should generate index values that are reasonably distributed throughout the range of possible values (such as in a SND)”

David's question was whether there is a way of showing that interestingness measures should be normally distributed. In fact, this question raises a number of important issues, and this post explores a few of them.

First, as I noted in my comment in response to this question, the short answer is that the “standard behavior” we expect from most estimators – i.e., most characterizations computed from a set of uncertain data values – is asymptotic normality. That is, for sufficiently large data samples, we expect most of the standard characterizations we might compute – variance estimates, correlation estimates, least squares regression coefficients, or linear model parameters computed using robust estimators like Least Trimmed Squares (LTS) – to be approximately normally distributed. I discuss this idea in Chapter 6 of Exploring Data in Engineering, the Sciences, and Medicine, beginning with the Central Limit Theorem (CLT), which says, very roughly, that “averages of N numbers tend to Gaussian limits as N becomes infinitely large.” There are exceptions – cases where asymptotic normality does not hold – but the phenomenon is common enough to make violations noteworthy. In fact, one of the reasons the conceptually simpler least median of squares (LMS) estimator was dropped in favor of the more complex least trimmed squares (LTS) estimator is that the LMS estimator is not asymptotically normal, approaching a non-normal limiting distribution at an anomalously slow rate. In contrast, the LTS estimator is asymptotically normal, with the standard convergence rate (i.e., the standard deviation of the estimator approaches zero inversely with the square root of the sample size as the sample becomes infinitely large).

The longer – and far less satisfying – answer to the question of how interestingness measures should be distributed is, “it depends,” as the following discussion illustrates.

One of the many cases where asymptotic normality has been shown to hold is Gini’s mean difference, which forms the basis for one of the four interestingness measures I discussed in my earlier post. While this result appears to provide a reasonably specific answer to the question of how this interestingness measure is distributed, there are two closely related issues that greatly complicate the matter. The first is a question that plagues many applications of asymptotic normality, which is quite commonly used for constructing approximate confidence intervals: how large must the sample be for this asymptotic approximation to be reasonable? (The confidence intervals described in my first post on odds ratios, for example, were derived on the basis of asymptotic normality.) Once again, the answer is, “it depends,” as the examples discussed here demonstrate. The second issue is that the interestingness measures that I describe in Exploring Data and that I discussed in my previous post are normalized measures. In particular, the Gini interestingness measure is defined as Gini’s mean difference divided by its maximum possible value, giving a measure that is bounded between 0 and 1. An important aspect of the Gaussian distribution is that it assigns a nonzero probability to any real value, although the probability associated with values more than a few standard deviations from the mean rapidly becomes very small. Thus, for a bounded quantity like the Gini interestingness measure, a Gaussian approximation can only be reasonable if the standard deviation is small enough that the probability of exhibiting values less than 0 or greater than 1 is acceptably small. Further, this “maximum reasonable standard deviation” will depend on the mean value of the estimator: if the computed interestingness measure is approximately 0.5, the maximum feasible standard deviation for a reasonable Gaussian approximation is considerably larger than if the computed interestingness measure is 0.01 or 0.99. As a practical matter, this usually means that, for a fixed sample size, the shape of the empirical distribution of normalized interestingness values for sequences exhibiting a given degree of heterogeneity (i.e., a specified “true” interestingness value) will vary significantly in shape (specifically, symmetry) as this “true” interestingness value varies from near 0 to approximately 0.5, to near 1. This idea is easily demonstrated in R on the basis of simulated examples.

The essential idea behind this simulation study is the following. Suppose we have a categorical variable that can assume any one of L distinct levels, and suppose we want to generate random samples of this variable where each level i has a certain probability p_i of occurring in our sample. The multinomial distribution discussed in Sec. 10.2.1 of Exploring Data is appropriate here, characterizing the number of times n_i we observe each of these levels in a random sample of size N, given the probability of observing each level i. The simulation strategy considered here then proceeds by specifying the number of levels (here, I take L = 4) and the probability of observing each level. Then, I generate a large number B of random samples (here, I take B = 1000), each generated from the appropriate multinomial distribution. In what follows, I consider the following four cases:

Case A: four equally-represented levels, each with probability 0.25;
Case B: one dominant level (p₁ = 0.97) and three rare levels (p₂ = p₃ = p₄ = 0.01);
Case C: three equally-represented majority levels (p₁ = p₂ = p₃ = 0.33) and one rare minority level (p₄ = 0.01);
Case D: four distinct probabilities (p₁ = 0.05, p₂ = 0.10, p₃ = 0.25, p₄ = 0.60).

To generate these sequences, I use R’s built-in multinomial random number generator, rmultinom, with parameters n = 1000 (this corresponds to B in the above discussion: I want to generate 1000 multinomial random vectors), size = 100 (initially, I take each vector to be of length N = 100; later, I consider larger vectors), and, for Case A, I specify prob = c(0.25, 0.25, 0.25, 0.25), corresponding to the four p_i values specified above (for the other three cases, I specify the prob parameter appropriately). The value returned by this procedure is a matrix, with one row for each level and one column for each simulation; for each column, the number in the i^th row is n_i, the number of times the i^th level is observed in that simulated response. Dividing these counts by the sample size gives the fractional representation of each level, from which the normalized interestingness measures are then computed. For the Gini interestingness measure and Case A, these steps are accomplished using the following R code:

set.seed(101)

pvectA <- c(0.25, 0.25, 0.25, 0.25)

CountsA <- rmultinom(n = 1000, size = 100, prob = pvectA)

FractionsA <- CountsA/100

GiniA <- apply(FractionsA, MARGIN=2, gini.proc)

The first line of this sequence sets the seed value to initialize the random number generator: if you don’t do this, running the same procedure again will give you results that are statistically similar but not exactly the same as before; specifying the seed guarantees you get exactly the same results the next time you execute the command sequence. The second line defines the “true” heterogeneity of each sequence (i.e., the underlying probability of observing each level that the random number generator uses to simulate the data). The third line invokes the multinomial random number generator to return the desired matrix of counts, where each column represents an independent simulation and the rows correspond to the four possible levels of the response. The fourth line converts these counts to fractions by dividing by the size of each sample, and the last line uses the R function apply to invoke the Gini interestingness procedure gini.proc for each column of the fractions matrix: setting MARGIN = 2 tells R to “apply the indicated function gini.proc to the columns of FractionsA.” (It would also be possible to write this procedure as a loop in R, but the apply function is much faster.) Replacing gini.proc by bray.proc in the above sequence of commands would generate the Bray interestingness measures, using shannon.proc would generate the Shannon interestingness measures, and using simpson.proc would generate the Simpson interestingness measures (note that all of these procedures are available from the companion website for Exploring Data).

For this first example – i.e., Case A, where all four levels are equally represented – note that the “true” value for any of the four normalized interestingness measures I discussed in my earlier post is 0, the smallest possible value. Since the random samples generated by the multinomial random number generator rarely have identical counts for the four levels, the interestingness measure computed from each sample is generally larger than 0, but it can never be smaller than this value. This observation suggests two things. The first is that the average of these interestingness values is larger than 0, implying that any of the normalized measures I considered exhibit a positive bias for this case. The second consequence of this observation is that the distribution of any of these interestingness values is probably asymmetric. The first of these conclusions is evident in the plot below, which shows nonparametric density estimates for the Gini measure computed from the 1000 random simulations generated for each case. The left-most peak corresponds to Case A and it is clear from this plot that the average value is substantially larger than zero (specifically, the mean of these simulations is 0.114). This plot also suggests asymmetry, although the visual evidence is not overwhelming.

The right-most peak in the above plot corresponds to Case B, for which the “true” Gini interestingness value is 0.96. Here, the peak is much narrower, although again the average of the simulation values (0.971) is slightly larger than the true value, which is represented by the right-most vertical dashed line. Case C is modeled on the CapSurf variable from the UCI mushroom dataset that I discussed in my earlier post on interestingness measures: the first three levels are equally distributed, but there is also an extremely rare fourth level. The “true” Gini measure for this case is 0.32, corresponding to the second-from-left dashed vertical line in the above plot, and – as with the previous two cases – the distribution of Gini values computed from the multinomial random samples is biased above this correct value. Here, however, the shape of the distribution is more symmetric, suggesting possible conformance with the Gaussian expectations suggested above. Finally, Case D is the third peak from the left in the above plot, and here the distribution appears reasonably symmetric around the correct value of 0.6, again indicated by a vertical dashed line.

The above plot shows the corresponding results for the Shannon interestingness measure for the same four cases. Here, the distributions for the different cases exhibit a much wider range of variation than they did for the Gini measure. As before, Case D (the smoothest of the four peaks) appears reasonably symmetric in its distribution around the “true” value of 0.255, possibly consistent with a normal distribution, but this is clearly not true for any of the other three cases. The asymmetry seen for Case A (the left-most peak) appears fairly pronounced, the density for Case C (the second-from-left peak) appears to be multi-modal, and the right-most peak (Case B) exhibits both asymmetry and multimodal character. The key point here is that the observed distribution of interestingness values depends on both the specific simulation case considered (i.e., the true probability distribution of the levels), and the particular interestingness measure chosen.

As noted above, if we expect approximate normality for a data characterization on the basis of its asymptotic behavior, an important practical question is whether we are working with large enough samples for asymptotic approximations to be at all reasonable. Here, the sample size is N = 100: we are generating and characterizing a lot more sequences than this, but each characterization is based on 100 data observations. In the UCI mushroom example I discussed in my previous post, the sample size was much larger: that dataset characterizes 8,124 different mushrooms, so the interestingness measures considered there were based on a sample size of N = 8124. In fact, this can make an important difference, as the following example demonstrates. The above plot compares the Gini measure density estimates obtained from 1000 random samples, each of length N = 100 (the dashed curve), with those obtained from 1000 random samples, each of length N = 8124 (the solid curve), for Case D. The vertical dotted line corresponds to the true Gini value of 0.6 for this case. Note that both of these densities appear to be symmetrically distributed around this correct value, but the distribution of values is much narrower when they are computed from the larger samples.

These results are consistent with our expectations of asymptotic normality for the Gini measure, at least for cases that are not too near either extreme. Further confirmation of these general expectations for a “well-behaved interestingness measure” is provided by the four Q-Q plot shown below, generated using the qqPlot command in the car package in R that I discussed in a previous post. These plots compare the Shannon measures computed from Cases C (upper two plots) and D (lower two plots), for the smaller sample size N = 100 (left-hand plots) and the larger sample size N = 8124 (right-hand plots). In particular, in the upper left Q-Q plot (Case C, N = 100), compelling evidence of non-normality is given by both the large number of points falling outside the 95% confidence intervals for this Q-Q plot, and the “kinks” in the plot that reflect the multimodality seen in the second peak in the second plot above. In contrast, all of the points in the upper right Q-Q plot (Case C, N = 8124) fall within the 95% confidence limits for this plot, suggesting that the approximate normality assumption is reasonable for the larger sample. Qualitatively similar behavior is seen for Case D in the two lower plots: the results computed for N = 100 exhibit some evidence for asymmetry, with both the lower and upper tails lying above the upper 95% confidence limit for the Q-Q plot, while the results computed for N = 8124 show no such deviations.

It is important to emphasize, however, that approximate normality is not always to be expected. The plot below shows Q-Q plots for the Gini measure (upper plots) and the Shannon measure (lower plots) for Case A, computed from both the smaller samples (N = 100, left-hand plots) and the larger samples (N = 8124, right-hand plots). Here, all of these plots provide strong evidence for distributional asymmetry, consistent with the point made above that, while the computed interestingness value can exceed the true value of 0 for this case, it can never fall below this value. Thus, we expect both a positive bias and a skewed distribution of values, even for large sample sizes, and we see both of these features here.

To conclude, then, I reiterate what I said at the beginning: the answer to the question of what distribution we should expect for a “good” interestingness measure is, “it depends.” But now I can be more specific. The answer depends, first, on the sample size from which we compute the interestingness measure: the larger the sample, the more likely this distribution is to be approximately normal. This is, essentially, a restatement of the definition of asymptotic normality, but it raises the key practical question of how large N must be for this approximation to be reasonable. That answer also depends, on at least two things: first, the specific interestingness measure we consider, and second, the true heterogeneity of the data sample from which we compute it. Both of these points were illustrated in the first two plots shown above: the shape of the estimated densities varied significantly both with the case considered (i.e., the “true” heterogeneity) and whether we considered the Gini measure or the Shannon measure. Analogous comments apply to the other interestingness measures I discussed in my previous post.

Computing Odds Ratios in R

2011-05-07T08:21:00.000-07:00

In my last post, I discussed the use of odds ratios to characterize the association between edibility and binary mushroom characteristics for the mushrooms characterized in the UCI mushroom dataset. I did not, however, describe those computations in detail, and the purpose of this post is to give a brief discussion of how they were done. That said, I should emphasize that the primary focus of this blog is not R programming per se, but rather the almost limitless number of things that R allows you to do in the realm of exploratory data analysis. For more detailed discussions of the mechanics of R programming, two excellent resources are The R Book by Michael J. Crawley, and Modern Applied Statistics with S (Statistics and Computing) by W.N. Venables and B.D. Ripley.

As a reminder, I discuss the odds ratio in Chapter 13 of Exploring Data in Engineering, the Sciences, and Medicine, which may be viewed as an association measure between binary variables. As I discussed last time, for two binary variables, x and y, each taking the values 0 and 1, the odds ratio is defined on the basis of the following four numbers:

N₀₀ = the number of data records with x = 0 and y = 0

N₀₁ = the number of data records with x = 0 and y = 1

N₁₀ = the number of data records with x = 1 and y = 0

N₁₁ = the number of data records with x = 1 and y = 1

Specifically, the odds ratio is given by the following expression:

OR = N₀₀ N₁₁ / N₀₁ N₁₀

Similarly, confidence intervals for the odds ratio are easily constructed by appealing to the asymptotic normality of log OR, which has a limiting variance given by the square root of the sum of the reciprocals of these four numbers. The R procedure oddsratioWald.proc available from the companion website for Exploring Data computes the odds ratio and the upper and lower confidence limits at a specified level alpha from these four values:

oddsratioWald.proc <- function(n00, n01, n10, n11, alpha = 0.05){

# Compute the odds ratio between two binary variables, x and y,

# as defined by the four numbers nij:

# n00 = number of cases where x = 0 and y = 0

# n01 = number of cases where x = 0 and y = 1

# n10 = number of cases where x = 1 and y = 0

# n11 = number of cases where x = 1 and y = 1

OR <- (n00 * n11)/(n01 * n10)

# Compute the Wald confidence intervals:

siglog <- sqrt((1/n00) + (1/n01) + (1/n10) + (1/n11))

zalph <- qnorm(1 - alpha/2)

logOR <- log(OR)

loglo <- logOR - zalph * siglog

loghi <- logOR + zalph * siglog

ORlo <- exp(loglo)

ORhi <- exp(loghi)

oframe <- data.frame(LowerCI = ORlo, OR = OR, UpperCI = ORhi, alpha = alpha)

oframe

}

Including “alpha = 0.05” in the parameter list fixes the default value for alpha at 0.05, which yields the 95% confidence intervals for the computed odds ratio, based on the Wald approximation described above. An important practical point is that these intervals become infinitely wide if any of the four numbers N_ij are equal to zero; also, note that in this case, the computed odds ratio is either zero or infinite. Finally, it is worth noting that if the numbers N_ij are large enough, the procedure just described can encounter numerical overflow problems (i.e., the products in either the numerator or the denominator become too large to be represented in machine arithmetic). If this is a possibility, a better alternative is to regroup the computations as follows:

OR = (N₀₀ / N₀₁) x (N₁₁ / N₁₀)

To use the routine just described, it is necessary to have the four numbers defined above, which form the basis for a two-by-two contingency table. Because contingency tables are widely used in characterizing categorical data, these numbers are easily computed in R using the table command. As a simple example, the following code reads the UCI mushroom dataset and generates the two-by-two contingency table for the EorP and GillSize attributes:

> mushrooms <- read.csv("mushroom.csv")

> table(mushrooms$EorP, mushrooms$GillSize)

b n

e 3920 288

p 1692 2224

(Note that the first line reads the csv file containing the mushroom data; for this command to work as shown, it is necessary for this file to be in the working directory. Alternatively, you can change the working directory using the setwd command.)

To facilitate the computation of odds ratios, the following preliminary procedure combines the table command with the oddsratioWald.proc procedure, allowing you to compute the odds ratio and its level-alpha confidence interval from the two-level variables directly:

TableOR.proc00 <- function(x,y,alpha=0.05){

xtab <- table(x,y)

n00 <- xtab[1,1]

n01 <- xtab[1,2]

n10 <- xtab[2,1]

n11 <- xtab[2,2]

oddsratioWald.proc(n00,n01,n10,n11,alpha)

}

The primary disadvantage of this procedure is that it doesn’t tell you which levels of the two variables are being characterized by the computed odds ratio. In fact, this characterization describes the first level of each of these variables, and the following slight modification makes this fact explicit:

TableOR.proc <- function(x,y,alpha=0.05){

xtab <- table(x,y)

n00 <- xtab[1,1]

n01 <- xtab[1,2]

n10 <- xtab[2,1]

n11 <- xtab[2,2]

outList <- vector("list",2)

outList[[1]] <- paste("Odds ratio between the level [",dimnames(xtab)[[1]][1],"] of the first variable and the level [",dimnames(xtab)[[2]][1],"] of the second variable:",sep=" ")

outList[[2]] <- oddsratioWald.proc(n00,n01,n10,n11,alpha)

outList

}

Specifically, I have used the fact that the dimension names of the 2x2 table xtab correspond to the levels of the variables x and y, and I have used the paste command to include these values in a text string displayed to the user. (I have enclosed the levels in square brackets to make them stand out from the surrounding text, particularly useful here since the levels are coded as single letters.) Applying this procedure to the mushroom characteristics EorP and GillSize yields the following results:

> TableOR.proc(mushrooms$EorP, mushrooms$GillSize)

[[1]]

[1] "Odds ratio between the level [ e ] of the first variable and the level [ b ] of the second variable:"

[[2]]

LowerCI OR UpperCI alpha

1 15.62615 17.89073 20.48349 0.05

Almost certainly, the formatting I have used here could be improved – probably a lot – but the key point is to provide a result that is reasonably complete and easy to interpret.

Finally, I noted in my last post that if we are interested in using odds ratios to compare or rank associations, it is useful to code the levels so that the computed odds ratio is larger than 1. In particular, note that applying the above procedure to characterize the relationship between edibility and the Bruises characteristic yields:

> TableOR.proc(mushrooms$EorP, mushrooms$Bruises)

[[1]]

[1] "Odds ratio between the level [ e ] of the first variable and the level [ f ] of the second variable:"

[[2]]

LowerCI OR UpperCI alpha

1 0.09014769 0.1002854 0.1115632 0.05

It is clear from these results that both Bruises and GillSize exhibit odds ratios with respect to mushroom edibility that are significantly different from the neutral value 1 (i.e., the 95% confidence interval excludes the value 1 in both cases), but it is not obvious which variable has the stronger association, based on the available data. The following procedure automatically restructures the computation so that the computed odds ratio is larger than or equal to 1, allowing us to make this comparison:

AutomaticOR.proc <- function(x,y,alpha=0.05){

xtab <- table(x,y)

n00 <- xtab[1,1]

n01 <- xtab[1,2]

n10 <- xtab[2,1]

n11 <- xtab[2,2]

rawOR <- (n00*n11)/(n01*n10)

if (rawOR < 1){

n01 <- xtab[1,1]

n00 <- xtab[1,2]

n11 <- xtab[2,1]

n10 <- xtab[2,2]

iLevel <- 2

}

else{

iLevel <- 1

}

outList <- vector("list",2)

outList[[1]] <- paste("Odds ratio between the level [",dimnames(xtab)[[1]][1],"] of the first variable and the level [",dimnames(xtab)[[2]][iLevel],"] of the second variable:",sep=" ")

outList[[2]] <- oddsratioWald.proc(n00,n01,n10,n11,alpha)

outList

}

Note that this procedure first constructs the 2x2 table on which everything is based and then computes the odds ratio in the default coding: if this value is smaller than 1, the coding of the second variable (y) is reversed. The odds ratio and its confidence interval are then computed and the levels of the variables used in computing it are presented as before. Applying this procedure to the Bruises characteristic yields the following result, from which we can see that GillSize appears to have the stronger association, as noted last time:

> AutomaticOR.proc(mushrooms$EorP, mushrooms$Bruises)

[[1]]

[1] "Odds ratio between the level [ e ] of the first variable and the level [ t ] of the second variable:"

[[2]]

LowerCI OR UpperCI alpha

1 8.963532 9.971541 11.09291 0.05

Well, that’s it for now. Next time, I will come back to the question of how we should expect interestingness measures to be distributed and whether those expectations are met in practice.

Measuring association using odds ratios

2011-04-23T11:13:00.000-07:00

In my last two posts, I have used the UCI mushroom dataset to illustrate two things. The first was the use of interestingness measures to characterize categorical variables, and the second was the use of binary confidence intervals to visualize the relationship between a categorical predictor variable and a binary response variable. This second approach can be applied to categorical predictors having any number of levels, but in the case of a binary (i.e., two-level) predictor, an attractive alternative is to measure their association with odds ratios. The objective of this post is to illustrate this idea and highlight a few important details.

The above plots show the binomial confidence intervals discussed last time for four different binary mushroom characteristics: GillSize (upper left), GillAtt (upper right), Bruises (lower left), and StalkShape (lower right). Specifically, these plots show the estimated probability that mushrooms with each of the two possible values for these variables are edible. Thus, the upper left plot shows that mushrooms with GillSize characteristic “b” (“broad”) are much more likely to be edible than mushrooms with GillSize characteristic “n” (“narrow”). The other three plots have analogous interpretations: mushrooms with GillAtt value “a” (“attached”) are more likely to be edible than those with value “f” (“free”), mushrooms with bruises (Bruises value “t”) are more likely to be edible than those without (Bruises value “f”), and mushrooms with StalkShape value “t” (“tapering”) are slightly more likely to be edible than those with value “e” (“enlarging”). Also, while the smaller slopes for GillAtt and StalkShape suggest this association is weaker for these variables than for GillSize, where the slope appears much larger, it would be nice to have a quantitative measure of this degree of association that we could compare directly. This is particularly the case for GillSize and Bruises, where both associations appear to be reasonably strong, but since the reference lines run in opposite directions on the plots, it is difficult to reliably compare the slopes on the basis of appearance alone.

The odds ratio provides a simple quantitative association measure for these variables that allows us to make these comparisons directly. I discuss the odds ratio in Chapter 13 of Exploring Data in Engineering, the Sciences, and Medicine in connection with the practical implications of data type (e.g., numerical versus categorical data). The odds ratio may be viewed as an association measure between binary variables, and it is defined as follows. For simplicity, suppose x and y are two binary variables of interest and assume that they are coded so that they each take the values 0 or 1 – this assumption is easily relaxed, as discussed below, but it simplifies the basic description of the odds ratio. Next, define the following four numbers:

N₀₀ = the number of data records with x = 0 and y = 0

N₀₁ = the number of data records with x = 0 and y = 1

N₁₀ = the number of data records with x = 1 and y = 0

N₁₁ = the number of data records with x = 1 and y = 1

The odds ratio is defined in terms of these four numbers as

OR = N₀₀ N₁₁ / N₀₁ N₁₀

Since all of the four numbers appearing in this ratio are nonnegative, it follows that the odds ratio is also nonnegative and can assume any value between 0 and positive infinity. Further, if x and y are two statistically independent binary random variables, it can be shown that the odds ratio is equal to 1. Values greater than 1 imply that records with y = 1 are more likely to have x = 1 than x = 0, and similarly, that records with y = 0 are more likely to have x = 0 than x = 1; in other words, OR > 1 implies that the variables x and y are more likely to agree than they are to disagree. Conversely, odds ratio values less than 1 imply that the variables x and y are more likely to disagree: records with y = 1 are more likely to have x = 0 than x = 1, and those with y = 0 are more likely to have x = 1 than x = 0.

Often – as in the mushroom dataset – the binary variables are not coded as 0 or 1, but instead as two different categorical values. As a specific example, the binary response variable considered last time – the edibility variable EorP – assumes the values “e” (for “edible”) or “p” (for “poisonous” or “non-edible”). In the results presented here, we recode EorP to have the values 1 for edible mushrooms and 0 for non-edible mushrooms. For the mushroom characteristic GillSize shown in the upper left plot above, suppose we initially code the value “b” (“broad”) as 0 and the value “n” (“narrow”) as 1. This choice is arbitrary – we could equally well code “b” as 1 and “n” as zero – and its practical consequences are explored further below. For the coding just described, the odds ratio between mushroom edibility (EorP) and gill size (GillSize) is 0.056. Since this number is substantially smaller than 1, it suggests that edible mushrooms (y = 1) are unlikely to be associated with narrow gills (x = 1), a result that is consistent with the appearance of the upper left plot above.

An important practical issue in interpreting odds ratios is that of how much smaller or larger than 1 the computed odds ratio should be to be regarded as evidence for a “significant” association between the variables x and y. That is, since we are computing this ratio from uncertain data, we need a measure of precision for the odds ratio, like the binomial confidence intervals discussed in my last post: e.g., how much does the odds ratio change if some mushrooms previously declared edible is reclassified as poisonous, or if some additional mushrooms are added to our dataset? Fortunately, confidence intervals for the odds ratio are easily constructed. In his book Categorical Data Analysis (Wiley Series in Probability and Statistics), Alan Agresti notes that confidence intervals for the odds ratio can be computed directly by appealing to the fact that the odds ratio estimator is asymptotically normal, approaching a Gaussian distribution in the limit of large sample sizes. He does not give explicit results for these direct confidence intervals, however, because he does not recommend them. Instead, Agresti advocates the construction of confidence intervals for the log of the odds ratio and transforming them back to get upper and lower confidence limits for the odds ratio itself. This recommendation rests primarily on three practical points: first, that the log of the odds ratio approaches normality faster than the odds ratio itself does, so this approach yields more accurate confidence intervals; second, this approach guarantees a positive lower confidence limit for the odds ratio, which is not the case for the direct approach; and, third, the same result can be used to compute confidence intervals for both the odds ratio and its reciprocal, a result that is again not true for the direct approach and that will be useful in the discussion presented below.

For the gill size example, Agresti’s recommended procedure yields a 95% confidence interval between 0.049 and 0.064. Since this interval does not include the value 1, we conclude that there is evidence to support an association between a mushroom’s gill size and its edibility, at least for mushrooms in the UCI dataset. Applying this procedure to the GillAtt characteristic shown in the upper right plot above yields an estimated odds ratio of 0.097 with a 95% confidence interval between 0.059 and 0.157. Again, the fact that this interval does not include 1 supports the idea that the GillAtt characteristic is associated with edibility (again, for the mushrooms considered here), but the fact that this odds ratio is larger (i.e., closer to the neutral value 1) also suggests that this association is weaker than that between edibility and the GillShape characteristic. Again, this result is in agreement with the visual appearance of the upper right plot above, relative to that of the upper left plot. The advantage of the odds ratio over these plots is that it provides a quantitative measure that can be used to make more objective comparisons, removing the subjective visual judgment required in comparing plots.

Applying this procedure to the Bruises variable yields an odds ratio of 9.972, with a 95% confidence interval from 8.963 to 11.093. The fact that these values are larger than 1 implies that mushrooms whose bruise characteristics have been coded as 1 (here, “t” for “true” or “bruised”) are more likely to be edible than those whose characteristics have been coded as 0 (here, “f” for “false” or “not bruised”). As noted above, this coding is arbitrary, as were the earlier assignments. An extremely useful observation is that if we reverse this assignment – i.e., for this example, if we code “bruised” as 0 and “not bruised” as 1 – we simply exchange the numbers N₀₀ with N₁₀ and also the numbers N₁₁ with N₀₁. The effect of these exchanges on the odds ratio is a reciprocal transformation:

OR = N₀₀N₁₁/N₀₁N₁₀ -> N₁₀N₀₁/N₁₁N₀₀ = 1/OR

This observation provides a simple basis for comparing results like those for GillSize where the odds ratio is less than 1 with those for Bruises where the odds ratio is greater than 1. As with the visual comparisons discussed above, it is not obvious from the odds ratios computed so far which of these variables is more strongly associated with mushroom edibility. Reversing the coding for GillSize so that “b” is coded as 1 and “n” is coded as 0 changes the odds ratio from 0.059 to 1/0.059 = 17.857. Since this number is larger than the odds ratio of 9.972 for Bruises, we can conclude that GillSize is more strongly associated with edibility – i.e., it is a better predictor of edibility for the mushrooms considered here – than Bruises, at least for this dataset.

In fact, the same trick can be applied to the confidence intervals, illustrating the third advantage noted above for Agresti’s preferred approach to constructing these intervals. Specifically, the asymptotically normal approximation says that the log of the odds ratio has a mean of log OR and a standard deviation S that can be simply computed from the four numbers N₀₀, N₀₁, N₁₀, and N₁₁. Since the Gaussian distribution is symmetric about its mean and log(1/OR) = - log(OR), it follows that the log of the reciprocal odds ratio has the same approximate standard deviation S as log(OR). In practical terms, this means that if we reverse the coding of our binary predictor variables, it is a simple matter to compute new confidence intervals as follows:

New lower CI = 1/Old upper CI

New odds ratio = 1/Old odds ratio

New upper CI = 1/Old lower CI

(Note here that because the reciprocal transformation is order-reversing, the transformation of the lower confidence limit yields the new upper confidence limit, and vice-versa; for a more detailed discussion of order-preserving and order-reversing transformations in general and the reciprocal transformation in particular, refer to Chapter 12 of Exploring Data.)

Applying these transformation results to the odds ratio for Bruises yields a new odds ratio of 0.100, with a 95% confidence interval from 0.090 to 0.112, which we can now compare with the earlier results for GillSize and GillAtt. Alternatively, if we reverse the coding of the results for GillSize and GillAtt, we obtain odds ratios that are larger than 1 between edibility and “the more edible value” of each of these mushroom characteristics. This has the advantage of giving us a sequences of odds ratios, all larger than 1, with the largest value suggestive of the strongest association between each mushroom characteristic variable and edibility. For the four mushroom characteristics shown in the above four plots, this approach yields the following odds ratios and their 95% confidence intervals:

GillSize: Lower CI = 15.625, OR = 17.857, Upper CI = 20.408

GillAtt: Lower CI = 6.369, OR = 10.309, Upper CI = 16.949

Bruises: Lower CI = 8.963, OR = 9.972, Upper CI = 11.093

StalkShape: Lower CI = 1.384, OR = 1.512, Upper Ci = 1.651

These results suggest that, of these four variables, the best predictor of mushroom edibility is GillSize, followed by GillAtt as second-best, then Bruises, and finally StalkShape as least predictive. These conclusions are probably the same as we would draw based on a careful comparison of the plots shown above, but the odds ratios computed in the way just described lead us to these conclusions much more directly.

Finally, it is important to make three points. First, as I have noted before – but the point is important enough to bear repeating – the associations described here between these binary mushroom characteristics and edibility are based entirely on the UCI mushroom dataset. Thus, these conclusions are only as representative of mushrooms in general or in any particular setting as the UCI mushroom dataset is representative of this larger and/or different mushroom population. In particular, mushrooms from other locales or unusual environments may exhibit different relationships between edibility and gill size or other characteristics than the UCI mushrooms do. The second key point is that the results presented here only attempt to assess the predictability of a single binary mushroom characteristic in isolation. To get a more complete picture of the relationship between mushroom characteristics and edibility, it is necessary to explore more general multivariate analysis techniques like logistic regression. More about that later. Last but not least, the third point is that I realized in reviewing this post before I issued it that I hadn't included any actual R code to compute odds ratios. In my next post, I will remedy this problem, giving a detailed view of how the numbers presented here were obained.

Screening for predictive characteristics … and a mea culpa

2011-04-12T16:46:00.000-07:00

In my last post, I considered the UCI mushroom dataset and characterized the variables included there using four different interestingness measures. When I began drafting this post, my intention was to consider the question of how the different mushroom characteristics included in this dataset relate to each mushroom’s classification as edible or poisonous. In fact, I do consider this problem here, but in the process of working out the example, I discovered a minor typographical error in Exploring Data in Engineering, the Sciences, and Medicine that has somewhat less minor consequences. Specifically, in Eq. (9.67) on page 413, two square roots were omitted, making the result incorrect as stated. (More specifically, the term in curly brackets that appears twice should have an exponent of ½, like the different term that appears twice in curly brackets in Eq. (9.66) just above it on the same page.) The consequence of this omission is that the confidence intervals defined by Eq. (9.67) are too narrow; further, since this equation was used to implement the R procedure binomCI.proc available from the companion website, the results generated by this procedure are also incorrect. I have brought these errors to Oxford’s attention and have asked them to replace the original R procedure with a corrected update, but if you have already downloaded this procedure, you need to be aware of the missing square root. The rest of this post carries out my original plan – which was to show how binomial confidence intervals can be useful in screening categorical variables for their ability to predict a binary outcome like edibility.

Recall that the UCI mushroom dataset gives 23 characteristics for each of 8,124 mushrooms, including a binary classification of each mushroom as “edible” or “poisonous.” The question considered here is which – if any – of the 22 mushroom attributes included in the dataset is potentially useful in predicting edibility. The basic idea is the following: each of these predictors is a categorical variable that can take on any one of a fixed set of possible values, so we can examine the groups of mushrooms defined by each of these values and estimate the probability that the mushrooms in the group are edible. As a specific example, the mushroom characteristic CapSurf has four possible values: “f” (fibrous), “g” (grooves), “s” (smooth), or “y” (scaly). In this case, we want to estimate the probability that mushrooms with CapSurf = f are edible, the probability that those with CapSurf = g are edible, and similarly for CapSurf = s and y. The most common approach to this problem is to estimate these probabilities as the fractions of edible mushrooms in each group: P_edible = n_edible/n_group. The difficulty with this number, taken by itself, is that it doesn’t tell us how much weight to give the final result: if we have one edible mushroom in a group of five, we get P_edible = 0.200, and we get the same result if we have 200 edible mushrooms in a group of 1,000. We are likely to put more faith in the second result than in the first, however, because it has a lot more weight of evidence behind it. For example, if we add a single edible mushroom to each group, our probability estimate for the first case increases from P_edible = 0.200 to P_edible = 0.333, while in the second case, the estimated probability only increases to P_edible = 0.201. Even worse, if we remove one edible mushroom from the first group, P_edible drops to 0.000, while in the second case, it only drops to 0.199.

This is where statistics can come to our rescue: in addition to computing the point estimate P_edible of the probability that a mushroom is edible, we can also compute confidence intervals, which quantify the uncertainty in this result. That is, a confidence interval is defined as a set of values that has at least some specified probability of containing the true but unknown value of P_edible. A common choice is 95% confidence limits: the true value of P_edible lies between some lower limit P_- and some upper limit P₊ with probability at least 95%. One of the key points of this post is that these intervals can be computed in more than one way, and the way that was widely adopted as “the standard method” for a long time has been found to be inadequate. Fortunately, a simple alternative is available that gives much better results, at least if you implement it correctly. The details follow, all based on the material presented in Section 9.7 of Exploring Data.

This standard method relies on the assumption of asymptotic normality: for a “sufficiently large” group (i.e., for “large enough” values of n_group and possibly n_edible), the estimator P_edible should approaches a Gaussian limiting distribution with variance P_edible(1 – P_edible)/n_group. If we assume our sample is large enough for this to be a good approximation, we can rely on known results for the Gaussian distribution to construct our confidence intervals. As a specific example, the 95% confidence interval would be centered at P_edible with upper and lower limits lying approximately plus or minus 1.96 standard deviations from this value, where the standard deviation is just the square root of the variance given above. The plot below shows the results obtained by applying this strategy to the groups of mushrooms defined by the four possible values of the CapSurf variable. Specifically, the open circles in this plot correspond to the estimated probability P_edible that a mushroom from the group defined by each CapSurf value is edible. The downward-pointing triangles represent the upper 95% confidence limit for this value, and the upward-pointing triangles represent the lower 95% confidence limit for this value. The horizontal dotted line corresponds to the average fraction of edible mushrooms in the UCI dataset, giving us a frame of reference for assessing the edibility results for each individual CapSurf value. That is, points lying well above this average line represent groups of mushrooms that are more edible than average, while points lying well below this average line represent groups of mushrooms that are less edible than average. The result obtained for level “g” clearly illustrates one difficulty with this approach: this group is extremely small, containing only four mushrooms, none of which are classified as edible. Thus, not only is P_edible zero, its associated variance is also zero, giving us zero-width confidence intervals. In words, this result is suggesting that mushrooms with grooved cap surfaces are never edible and that we are quite certain of this, despite the fact that this conclusion is only based on four mushrooms. In contrast, we seem to be less certain about the probability that scaly (“y”) or smooth (“s”) mushrooms are edible, despite the fact that these results are based on groups of 3,244 and 2,556 mushrooms, respectively.

An alternative approach that gives more accurate confidence intervals and also overcomes this particular difficulty is one proposed by Brown, Cai, and DasGupta. The details are given in Exploring Data (with the exception of the error noted at the beginning of this post, I believe they are correct), and they are somewhat messy, so I won’t repeat them here, but the basic ideas are, first, to add positive offsets to both n_edible and n_group in computing the probability that a mushroom is edible, and second, to modify the expression for the variance. Both of these modifications depend explicitly on the confidence level considered (i.e., the offsets are different for 95% confidence intervals than they are for 99% confidence intervals, as are the variance modifications), and they become negligible in the limit as both n_edible and n_group become very large. To see the impact of these modifications, the plot below gives the modified 95% confidence intervals for the CapSurf data, in the same general format as before. Comparing this plot with the one above, it is clear that the most dramatic difference is that for level “g,” the grooved cap mushrooms: whereas the asymptotic result suggested the mushrooms in this group were poisonous with absolute certainty, the very wide confidence intervals for this group in the plot below reflect the fact that this result is only based on four mushrooms and, while none of these four are edible, the confidence intervals extend from essentially zero probability of being edible to almost the average probability for the complete dataset. Thus, while we can conclude from this plot that mushrooms with grooved cap surfaces appear less likely than average to be edible, the available evidence isn’t enough to make this argument too strongly. In contrast, the results for the mushrooms with scaly cap surfaces (“y”) or smooth cap surfaces (“s”) are essentially identical to those presented above, consistent with the much larger groups on which these results are based.

Before leaving this example, it is worth showing how the results are changed in light of my typographical error in Eq. (9.67) of Exploring Data. The missing square roots were omitted from terms that define the width of the confidence intervals, and these terms are numerically smaller than 1. Since if x is smaller than 1, the square root of x is larger than x but still smaller than 1, the effect of these omitted square roots is to make the resulting confidence intervals too narrow (i.e., we are using the value x rather than the larger square root of x that we should be using to determine the width of the confidence intervals). This error causes our results to appear more precise than they really are. This effect may be seen clearly in the plot below, which corresponds to the two plots discussed above, but with the confidence intervals based on the erroneous implementation of the estimator of Brown et al.

Finally, the figure below shows four plots, each in the same format as those discussed above, corresponding to the P_edible estimates obtained by applying the method of Brown et al. to four different mushroom characteristics. The upper left plot shows the results obtained for the mushroom characteristic “Odor,” which appears to be highly predictive of edibility. Careful examination of these results reveals that, for the mushrooms in the UCI dataset, those with odors characterized as “a” (almond) or “l” (anise) are always edible, those with odors characterized as “c” (creosote), “f” (foul), “m” (musty), “p” (pungent), “s” (spicy), or “y” (fishy) are always poisonous, and those with no odor are more likely to be edible than not, but they can still be poisonous. In contrast, CapShape (upper right plot) appears much less predictive: some values seem to be strongly associated with edibility (“b” or “s”), while the levels “f” and “x” seem to convey no information at all: the likelihood that these mushrooms are edible is essentially the same as that of the complete collection, without regard to CapShape. The lower left plot shows the corresponding results for StalkRoot, which suggest that levels “c,” “e,” and “r” are more likely to be edible than average, level “b” conveys no information, and mushrooms where StalkRoot values are missing are somewhat more likely to be poisonous (the class “?”). This result is somewhat distressing, raising the possibility that the missing values for this variable are not missing at random, but that there may be some systematic mechanism at work (e.g., is the StalkRoot characterization somehow more difficult for poisonous mushrooms?). Finally, the lower right plot shows the results for the binary characteristic GillSize: it appears that mushrooms with GillSize “n” (narrow) are much more likely to be poisonous than those with GillSize “b” (broad). Because both the response (i.e., edibility) and the candidate predictor GillSize are binary in this case, an alternative – and arguably better – approach to characterizing their relationship is in terms of odds ratios, which I will take up in my next post.

Interestingness Measures

2011-04-03T08:06:00.000-07:00

Probably because I first encountered them somewhat late in my professional life, I am fascinated by categorical data types. Without question, my favorite book on the subject is Alan Agresti’s Categorical Data Analysis (Wiley Series in Probability and Statistics), which provides a well-integrated, comprehensive treatment of the analysis of categorical variables. One of the topics that Agresti does not discuss, however, is that of interestingness measures, a useful quantitative characterization of categorical variables that comes from the computer science literature rather than the statistics literature. As I discuss in Chapter 3 of Exploring Data in Engineering, the Sciences, and Medicine, interestingness measures are essentially numerical characterizations of the extent to which a categorical variable is uniformly distributed over its range of possible values: variables where the levels are equally represented are deemed “less interesting” than those whose distribution varies widely across these levels. Many different interestingness measures have been proposed, and Hilderman and Hamilton give an excellent survey, describing 13 different measures in detail. (I have been unable to find a PDF version on-line, but the reference is R.J. Hilderman and H.J. Hamilton, “Evaluation of interestingness measures for ranking discovered knowledge,” in Proceedings of the 5^th Asia-Pacific Conference on Knowledge Discovery and Data Mining, D. Chueng, G.J. Williams, and Q. Li, eds., Hong Kong, April, 2001, pages 247 to 259.) In addition, the authors present five behavioral axioms for characterizing interestingness measures. In Exploring Data, I consider four normalized interestingness measures that satisfy the following three of Hilderman and Hamilton’s axioms:

1. The minimum value principle: the measure exhibits its minimum value when the variable is uniformly distributed over its range of possible values. For the normalized measures considered here, this means the measure takes the value 0.

2. The maximum value principle: the measure exhibits its maximum value when the variable is completely concentrated on one of its several possible values. For the normalized measures considered here, this means the measure takes the value 1.

3. The permutation-invariance principle: re-labeling the levels of the variable does not change the value of the interestingness measure.

To compute the interestingness measures considered here, it is necessary to first compute the empirical probabilities that a variable assumes each of its possible values. To be specific, assume that the variable x can take any one of M possible values, and let p_idenote the fraction of the N observed samples of x that assume the i^th possible value. All of the interestingness measures considered here attempt to characterize the extent to which these empirical probabilities are constant, i.e. the extent to which p_i is approximately equal to 1/M for all i. Probably the best known of the four interestingness measures I consider is the Shannon entropy from information theory, which is based on the average of the product p_i log p_iover all i. A second measure is the normalized version of Gini’s mean difference from statistics, which is the average distance that p_i lies from p_j for all i distinct from j, and a third – Simpson’s measure – is a normalized version of the variance of the p_i values. The fourth characterization considered in Exploring Data is Bray’s measure, which comes from ecology and is based on the average of the smaller of p_i and 1/M. The key point here is that, because these measures are computed in different ways, they are sensitive to different aspects of the distributional heterogeneity of a categorical variable over its range of possible values. Specifically, since all four of these measures assume the value 0 for uniformly distributed variables and 1 for variables completely concentrated on a single value, they can only differ for intermediate degrees of heterogeneity.

The R procedures bray.proc, gini.proc, shannon.proc, and simpson.proc are all available from the Exploring Data companion website, each implementing the corresponding interestingness measure. To illustrate the use of these procedures, they are applied here to the UCI Machine Learning Repository Mushroom dataset, which gives 23 categorical characterizations of 8,124 different species of mushrooms, taken from The Audubon Society Field Guide to North American Mushrooms (G.H. Lincoff, Pres., published by Alfred A. Knopf, New York, 1981). A typical characterization is gill color, which exhibits 12 distinct values, each corresponding to a one-character color code (e.g., “p” for pink, “u” for purple, etc.). To evaluate the four interestingness measures for any of these attributes, it is necessary to first compute its associated empirical probability vector. This is easily done with the following R function:

ComputePvalues <- function(x){

xtab <- table(x)

pvect <- xtab/sum(xtab)

pvect

}

Given this empirical probability vector, the four R procedures listed above can be used to compute the corresponding interestingness measure. As a specific example, the following sequence gives the values for the four interestingness measures for the seven-level variable “Habitat” from the mushroom dataset:

> x <- mushroom$Habitat

> pvect <- ComputePvector(x)

> bray.proc(pvect)

> 0.427212

> gini.proc(pvect)

> 0.548006

> shannon.proc(pvect)

> 0.189719

> simpson.proc(pvect)

> 0.129993

These results illustrate the point noted above, that while all of these measures span the same range – from 0 for completely uniform level distributions to 1 for variables completely concentrated on one level – in general, the values of the four measures are different. These differences arise from the fact that the different interestingness measures weight specific types of deviations from uniformity differently. To see this, note that the seven levels of the habitat variable considered here have the following distribution:

> table(mushroom$Habitat)

> d g l m p u w

> 3148 2148 832 292 1144 368 192

It is clear from looking at these numbers that the seven different habitat levels are not all equally represented in this dataset, with the most common level (“d”) occurring about 15 times as often as the rarest level (“w”). The average representation is approximately 1160, so the two most populous levels occur much more frequently than average, one level occurs with about average frequency (“p”), and the other four levels occur anywhere between half as often as average and one tenth as often. It is clear that the Gini measure is the most sensitive to these deviations from homogeneity, at least for this example, while the Simpson measure is the least sensitive. These observations raise the following question: to what extent is this behavior typical, and to what extent is it specific to this particular example? The rest of this post examines this question further.

The four plots shown above compare these different interestingness measures for all of the variables included in the UCI mushroom dataset. The upper left plot shows the Gini measure plotted against the Shannon measure for all 23 of these categorical variables. The dashed line in the plot is the equality reference line: if both measures were identical, all points would lie along this line. The fact that all points lie above this line means that – at least for this dataset – the Gini measure is always larger than the Shannon measure. The upper right plot shows the Simpson measure plotted against the Shannon measure, and here it is clear that the Simpson measure is generally smaller than the Shannon measure, but there are exceptions. Further, it is also clear from this plot that the differences between the Shannon and Simpson measures are substantially less than those between the Shannon and Gini measures, again, at least for this dataset. The lower left plot shows the Bray measure plotted against the Shannon measure, and this plot has the same general shape as the one above for the Gini versus Shannon measures. The differences between the Bray and Shannon measures appear slightly less than those between the Gini and Shannon measures, suggesting that the Bray measure is slightly smaller than the Gini measure. The lower right plot of the Bray versus the Gini measures shows that this is generally true: some of the points on this plot appear to fall exactly on the reference line, while the others fall slightly below this line, implying that the Bray measure is generally approximately equal to the Gini measure, but where they differ, the Bray measure is consistently smaller. Again, it is important to emphasize that the results presented here are based on a single dataset, and they may not hold in general, but these observations do what a good exploratory analysis should: they raise the question.

Another question raised by these results is whether there is anything special about the few points that violate the general rule that the Simpson measure is less than the Shannon measure. The above figure repeats the upper right plot from the four shown earlier, with the Simpson measure plotted against the Shannon measure. Here, the six points that correspond to binary characterizations are represented as solid circles, and they correspond to six of the eight points that fall above the line, implying that the Simpson measure is greater than the Shannon measure for these points. The fact that these are the only binary variables in the dataset suggests – but does not prove – that the Simpson measure tends to be greater than the Shannon measure for binary variables. Of the other two points above the reference line, one represents the only three-level characterization in the mushroom dataset (specifically, the point corresponding to “RingNum,” a categorical variable based on the number of rings on the mushroom, labeled on the plot). The last remaining point above the equality reference line corresponds to one of four four-level characterizations in the dataset; the other three four-level characterizations fall below the reference line, as do all of those with more than four levels. Taken together, these results suggest that the number of levels in a categorical variable influences the Shannon and Simpson measures differently. In particular, it appears that those cases where the Simpson measure exceeds the Shannon measure tend to be variables with relatively few levels.

In fact, these observations raise a subtle but extremely important point in dealing with categorical variables. The metadata describing the mushroom dataset indicates that the variable VeilType is binary, with levels universal (“u”) and partial (“p”), but the type listed for all 8,124 mushrooms in the dataset is “p.” As a consequence, VeilType appears in this analysis as a one-level variable. Since the basic numerical expressions for all four of the interestingness measures considered here all become indeterminate for one-level variables, this case is tested for explicitly in the R procedures used here, and it is assigned a value of 0: this seems reasonable, representing a classification of “fully homogeneous, evenly spread over the range of possible values,” as it is difficult to imagine declaring a one-level variable to be strongly heterogeneous or “interesting.” Here, this result is a consequence of defining the number of levels for this categorical variable from the data alone, neglecting the possibility of other values that could be present but are not. In particular, if we regard this variable as binary – in agreement with the metadata – all four of the interestingness measures would yield the value 1, corresponding to the fact that the observed values are fully concentrated in one of two possible levels. This example illustrates the difference between internal categorization – i.e., determination of the number of levels for a categorical variable from the observed data alone – and external categorization, where the number of levels is specified by the metadata. As this example illustrates, the difference can have important consequences for both numerical data characterizations and their interpretation.

Finally, the above figure shows the Gini measure plotted against the Shannon measure for all 23 of the mushroom characteristics, corresponding to the upper left plot in the first figure, but with three additions. First, the solid curve represents the result of applying the lowess nonparametric smoothing procedure to the scatterplot of Gini versus Shanon measures. Essentially, this curve fills in the line that your eye suggests is present, approximating the nonlinear relationship between these two measures that is satisfied by most of the points on the plot. Second, the solid circle represents the variable whose results lie farthest from this curve, which is CapSurf, a four-level categorical variable describing the surface of the mushroom cap as fibrous (“f”), grooved (“g”), scaly (“y”), or smooth (“s”), distributed as follows: 2,320 “f,” 4 “g,” 2,556 “y,” and 3224 “s.” Note that if the rare cagetory “g” were omitted from consideration, this variable would be almost uniformly distributed over the remaining three values, and this modified result (i.e., with the rare level "g" omitted from the analysis) corresponds to the solid square point that falls on the lowess curve at the lower left end. Thus, the original CapSurf result appears to represent a specific type of deviation from uniformity that violates the curvilinear relationship that seems to generally hold – at least approximately – between the Gini and Shannon interestingness measures.

The real point of these examples has been two-fold. First, I wanted to introduce the notion of an interestingness measure and illustrate the application of these numerical summaries to categorical variables, something that is not covered in typical introductory statistics or data analysis courses. Second, I wanted to show that – in common with most other summary statistics – different data characteristics influence alternative measures differently. Thus, while it appears that the Simpson measure is generally slightly smaller than the Shannon measure – at least for the examples considered here – this behavior depends strongly on the number of levels the categorical variable can assume, with binary variables consistently violating this rule of thumb (again, at least for the examples considered here). Similarly, while it appears that the Gini measure always exceeds the Shannon measure, often substantially so, there appear to be certain types of heterogeneous data distributions for which this difference is substantially smaller than normal. Thus, as is often the case, it can be extremely useful to compute similar characterizations of the same type and comapre them, since unexpected differences can illuminate features of the data that are not initially apparent and which may turn out to be interesting.

The Many Uses of Q-Q Plots

2011-03-23T18:17:00.000-07:00

My last four posts have dealt with boxplots and some useful variations on that theme. Just after I finished the series, Tal Galili, who maintains the R-bloggers website, pointed me to a variant I hadn’t seen before. It's called a beeswarm plot, and it's produced by the beeswarm package in R. I haven’t played with this package a lot yet, but it does appear to be useful for datasets that aren’t too large and that you want to examine across a moderate number of different segments. The plot shown below provides a typical illustration: it shows the beeswarm plot comparing the potassium content of different cereals, broken down by manufacturer, from the UScereal dataset included in the MASS package in R. I discussed this data example in my first couple of boxplot posts and I think this is a case where the beeswarm plot gives you a more useful picture of how the data points are distributed than the boxplots do. For more information about the beeswarm package, I recommend Tal's post. More generally, anyone interested in learning more about what you can do with the R software package should find the R-blogger website extremely useful.

Besides boxplots, one of the other useful graphical data characterizations I discuss in Exploring Data in Engineering, the Sciences, and Medicine is the quantile-quantile (Q-Q) plot. The most common form of this characterization is the normal Q-Q plot, which represents an informal graphical test of the hypothesis that a data sequence is normally distributed. That is, if the points on a normal Q-Q plot are reasonably well approximated by a straight line, the popular Gaussian data hypothesis is plausible, while marked deviations from linearity provide evidence against this hypothesis. The utility of normal Q-Q plots goes well beyond this informal hypothesis test, however, which is the main point of this post. In particular, the shape of a normal Q-Q plot can be extremely useful in highlighting distributional asymmetry, heavy tails, outliers, multi-modality, or other data anomalies. The specific objective of this post is to illustrate some of these ideas, expanding on the discussion presented in Exploring Data.

The above figure shows four different normal Q-Q plots that illustrate some of the different data characteristics these plots can emphasize. The upper left plot demonstrates that normal Q-Q plots can be extremely effective in highlighting glaring outliers in a data sequence. This plot shows the annual number of traffic deaths per ten thousand drivers over an unspecified time period, for 25 of the 50 states in the U.S., plus the District of Columbia. This plot was constructed from the road dataset included in the MASS package in R, which gives the numbers of deaths, the numbers of drivers (in tens of thousands), and several other characteristics for each of these regions. Based on the interpretation of normal Q-Q plots offered above, the normal distribution hypothesis appears fairly reasonable for this data sequence, in all cases except the point in the extreme upper right. This point corresponds to the state of Maine, which exhibited 26 deaths per ten thousand drivers, well above the average of approximately 5 for all other regions considered.

It is not clear why the reported traffic death rate is so high for Maine. The scatterplot above shows the reported traffic deaths for each state or district against the number of drivers, in tens of thousands. The dashed line in the plot corresponds to the average traffic death rate for all regions except Maine, and it is clear that this line fits most of the data points reasonably well, with Maine (the solid point) representing the most glaring exception. Although it still leaves us wanting to know more, this plot suggests that the number of deaths for Maine is unusually high, rather than the number of drivers being unusually low, which might be a more tempting explanation.

The Q-Q plot for this denominator variable – i.e., for the number of drivers – is shown as the upper right plot in the original set of four shown above. There, the fact that both tails of the distribution lie above the reference line is suggestive of distributional asymmetry, a point examined further below using Q-Q plots for other reference distributions. Also, note that both of the upper Q-Q plots shown above are based on only 26 data values, which is right at the lower limit on sample size that various authors have suggested for normal Q-Q plots to be useful (see the discussion of normal Q-Q plots in Section 6.3.3 of Exploring Data for details). The tricky issues of separating outliers, asymmetry, and other potentially interesting data characteristics in samples this small is greatly facilitated using the Q-Q plot confidence intervals discussed below.

The lower left Q-Q plot in the above sequence is that for the Old Faithful geyser dataset faithful included with the base R package. As I have discussed previously, the eruption duration data exhibits a pronounced bimodal distribution, which may be seen clearly in nonparametric density estimates computed from these data values. Normal Q-Q plots constructed from bimodal data typically exhibit a “kink” like the one seen in this plot. A crude way of explaining this behavior is the following: the lower portion of the Q-Q plot is very roughly linear, suggesting a very approximate Gaussian distribution, corresponding to the first mode of the eruption data distribution (i.e., the durations of the shorter group of eruptions). Similarly, the upper portion of the Q-Q plot is again very roughly linear, but with a much different intercept that corresponds to the larger mean of the second peak in the distribution (i.e., the durations of the longer group of eruptions). To connect these two “roughly linear” local segments, the curve must exhibit a “kink” or rapid transition region between them. By the same reasoning, more general multi-modal distributions will exhibit more than one such “kink” in their Q-Q plots. Finally, the lower right Q-Q plot in the collection above was constructed from the Pima Indians diabetes dataset available from the UCI Machine Learning Repository. This dataset includes a number of clinical measurements for 768 female members of the Pima tribe of Native Americans, including their diastolic blood pressure. The lower right Q-Q plot was constructed from this blood pressure data, and its most obvious feature is the prominent lower tail anomaly. In fact, careful examination of this plot reveals that these points correspond to the value zero, which is not realistic for any living person. What has happened here is that zero has been used to code missing values, both for this variable and several others in this dataset. This observation is important because the metadata associated with this dataset indicates that there is no missing data, and a number of studies in the classification literature have proceeded under the assumption that this is true. Unfortunately, this assumption can lead to badly biased results, a point discussed in detail in a paper I published in SIGKDD Explorations (Disguised Missing Data paper PDF). The point of the example presented here is to show that normal Q-Q plots can be extremely effective in highlighting this kind of data anomaly.

The normal Q-Q plots considered so far were constructed using the qqnorm procedure available in base R, and the reference lines shown in these plots were constructed using the qqline command. It is not difficult to construct Q-Q plots for other reference distributions using procedures in base R, but a much simpler alternative is to use the qqPlot command in the optional car package. This R add-on package was developed in association with the book An R Companion to Applied Regression, by Fox and Weisberg, and it includes a number of very useful procedures. The default options of the qqPot procedure automatically generate a reference line, along with upper and lower 95% confidence intervals for the plot, which are particularly useful for small samples like the road dataset. The figure below shows a normal Q-Q plot for the number of traffic deaths per 10,000 drivers generated using the qqPlot package. The fact that all of the points but the one obvious outlier fall within the 95% confidence limits suggest that the scatter around the reference line seen for these 25 observations is small enough to be consistent with a normal reference distribution. Further, these confidence limits also emphasize how much the outlying result for the state of Maine violates this normality assumption.

Another advantage of the qqPlot command is that it provides the basis for very easy generation of Q-Q plots for essentially any reference distribution that is available in R, including those available in add-on packages like gamlss.dist, which supports an extremely wide range of distributions (generalized inverse Gaussian distributions, anyone?). This capability is illustrated in the four Q-Q plots shown below, all generated with the qqPlot command for non-Gaussian distributions. In all of these plots, the data corresponds to the driver counts for the 26 states and districts summarized in the road dataset. Motivation for the specific Q-Q plots shown here is that the four distributions represented by these plots are all better suited to capturing the asymmetry seen in the normal Q-Q plot for this data sequence than the symmetric Gaussian distribution is. The upper left plot shows the results obtained for the exponential distribution which, like the Gaussian distribution, does not require the specification of a shape parameter. Comparing this plot with the normal Q-Q plot shown above for this data sequence, it is clear that the exponential distribution is more consistent with the driver data than the Gaussian distribution is. The data point in the extreme upper right does fall just barely outside the 95% confidence limits shown on this plot, and careful inspection reveals that the points in the lower left fall slightly below these confidence limits, which become quite narrow at this end of the plot.

The exponential distribution represents a special case of the gamma distribution, with a shape parameter equal to 1. In fact, the exponential distribution exhibits a J-shaped density, decaying from a maximum value at zero, and it corresponds to a “dividing line” within the gamma family: members with shape parameters larger than 1 exhibit unimodal densities with a single maximum at some positive value, while gamma distributions with shape parameters less than 1 are J-shaped like the exponential distribution. To construct Q-Q plots for general members of the gamma family, it is necessary to specify a particular value for this shape parameter, and the other three Q-Q plots shown above have done this using the qqPlot command. Comparing these plots, it appears that increasing the shape parameter causes the points in the upper tail to fall farther outside the 95% confidence limits, while decreasing the shape parameter better accommodates these upper tail points. Conversely, decreasing the shape parameter causes the cluster of points in the lower tail to fall farther outside the confidence limits. It is not obvious that any of the plots shown here suggest a better fit than the exponential distribution, but the point of this example was to show the flexibility of the qqPlot procedure in being able to pose the question and examine the results graphically. Alternatively, the Weibull distribution – which also includes the exponential distribution as a special case – might describe these data values better than any member of the gamma distribution family, and these plots can also be easily generated using the qqPlot command (just specify dist = “weibull” instead of dist = “gamma”, along with shape = a for some positive value of a other than 1).

Finally, one cautionary note is important here for those working with very large datasets. Q-Q plots are based on sorting data, something that can be done quite efficiently, but which can still take a very long time for a really huge dataset. As a consequence, while you can attempt to construct Q-Q plots for sequences of hundreds of thousands of points or more, you may have to wait a long time to get your plot. Further, it is often true that plots made up of a very large number of points reduce to ugly-looking dark blobs that can use up a lot of toner if you make the further mistake of trying to print them. So, if you are working with really enormous datasets, my suggestion is to construct Q-Q plots from a representative random sample of a few hundred or a few thousand points, not hundreds of thousands or millions of points. It will make your life a lot easier.