When are averages useless?

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

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

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. 

More about that in future posts.

A Brief Introduction to Mixture Distributions

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

  #

  y

}

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

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.