Implementing the CountSummary Procedure

In my last post, I described and demonstrated the CountSummary procedure to be included in the ExploringData package that I am in the process of developing.  This procedure generates a collection of graphical data summaries for a count data sequence, based on the distplot, Ord_plot, and Ord_estimate functions from the vcd package.  The distplot function generates both the Poissonness plot and the negative-binomialness plot discussed in Chapters 8 and 9 of Exploring Data in Engineering, the Sciences and Medicine.  These plots provide informal graphical assessments of the conformance of a count data sequence with the two most popular distribution models for count data, the Poisson distribution and the negative-binomial distribution.  As promised, this post describes the R code needed to implement the CountSummary procedure, based on these functions from the vcd package.

The key to this implementation lies in the use of the grid package, a set of low-level graphics primitives included in base R.  As I mentioned in my last post, the reason this was necessary - instead of using higher-level graphics packages like lattice or ggplot2 - was that the vcd package is based on grid graphics, making it incompatible with base graphics commands like those used to generate arrays of multiple plots.  The grid package was developed by Paul Murrell, who provides a lot of extremely useful information about both R graphics in general and grid graphics in particular on his home page, including the article “Drawing Diagrams with R,” which provides a nicely focused introduction to grid graphics.  The first example I present here is basically a composite of the first two examples presented in this paper.  Specifically, the code for this example is:

library(grid)

grid.newpage()

pushViewport(viewport(width = 0.8, height = 0.4))

grid.roundrect()

grid.text("This is text in a box")

popViewport()

The first line of this R code loads the grid package and the second tells this package to clear the plot window; failing to do this will cause this particular piece of code to overwrite whatever was there before, which usually isn’t what you want.  The third line creates a viewport, into which the plot will be placed.  In this particular example, we specify a width of 0.8, or 80% of the total plot window width, and a height of 0.4, corresponding to 40% of the total window height.  The next two lines draw a rectangular box with rounded corners and put “This is text in a box” in the center of this box.  The advantage of the grid package is that it provides us with simple graphics primitives to draw this kind of figure, without having to compute exact positions (e.g., in inches) for the different figure components.  Commands like grid.text provide useful defaults (i.e., put the text in the center of the viewport), which can be overridden by specifying positional parameters in a variety of ways (e.g., left- or right-justified, offsets in inches or lines of text, etc.).  The results obtained using these commands are shown in the figure below.

The code for the second example is a simple extension of the first one, essentially consisting of the added initial code required to create the desired two-by-two plot array, followed by four slightly modified copies of the above code.  Specifically, this code is:

grid.newpage()

pushViewport(viewport(layout=grid.layout(nrow=2,ncol=2)))

pushViewport(viewport(layout.pos.row=1,layout.pos.col=1))

grid.roundrect(width = 0.8, height=0.4)

grid.text("Plot 1 goes here")

popViewport()

pushViewport(viewport(layout.pos.row=1,layout.pos.col=2))

grid.roundrect(width = 0.8, height=0.4)

grid.text("Plot 2 goes here")

popViewport()

pushViewport(viewport(layout.pos.row=2,layout.pos.col=1))

grid.roundrect(width = 0.8, height=0.4)

grid.text("Plot 3 goes here")

popViewport()

pushViewport(viewport(layout.pos.row=2,layout.pos.col=2))

grid.roundrect(width = 0.8, height=0.4)

grid.text("Plot 4 goes here")

popViewport()

Here, note that the first “pushViewport” command creates the two-by-two plot array we want, by specifying “layout = grid.layout(nrow=2,ncol=2)”.  As in initializing a data frame in R, we can create an arbitrary two-dimensional array of grid graphics viewports – say m by n – by specifying “layout = grid.layout(nrow=m, ncol=n)”.  Once we have done this, we can use whatever grid commands – or grid-compatible commands, such as those generated by the vcd package – we want, to create the individual elements in our array of plots.  In this example, I have basically repeated the code from the first example to put text into rounded rectangular boxes in each position of the plot array.  The two most important details are, first, the “pushViewport” command at the beginning of each of these individual plot blocks specifies which of the four array elements the following plot will go in, and second, the “popViewport()” command at the end of each block, which tells the grid package that we are finished with this element of the array.  If we leave this command out, the next “pushViewport” command will not move to the desired plot element, but will simply overwrite the previous plot.  Executing this code yields the plot shown below.

The final example replaces the text in the above two-by-two example with the plots I want for the CountSummary procedure.  Before presenting this code, it is important to say something about the structure of the resulting plot and the vcd commands used to generate the different plot elements.  The first plot – in the upper left position of the array shown below – is an Ord plot, generated by the Ord_plot command, which does two things.  The first is to generate the desired plot, but the second is to return estimates of the intercept and slope of one of the two reference lines in the plot.  The first of these lines is fit to the points in the plot via ordinary least squares, while the second – the one whose parameters are returned – is fit via weighted least squares, to down-weight the widely scattered points seen in this plot that correspond to cases with very few observations.  The intent of the Ord plot is to help us decide which of several alternative distributions – including both the Poisson and the negative-binomial – fits our count data sequence better.   This guidance is based on the reference line parameters, and the Ord_estimate function in the vcd package transforms these parameter estimates into distributional recommendations and the distribution parameter values needed by the distplot function in the vcd package to generate either the Poissonness plot or the negative-binomialness plot for the count data sequence.  Although these recommendations are sometimes useful, it is important to emphasize the caution given in the vcd package documentation:

           

“Be careful with the conclusions from Ord_estimate as it implements just some simple heuristics!”

In the CountSummary procedure, I use these results both to generate part of the text summary in the upper right element of the plot array, and to decide which type of plot to display in the lower right element of this array.  Both this plot and the Poissonness reference plot in the lower left element of the display are created using the distplot command in the vcd package.  I include the Poissonness reference plot because the Poisson distribution is the most commonly assumed distribution for count data – analogous in many ways to the Gaussian distribution so often assumed for continuous-valued data – and, by not specifying the single parameter for this distribution, I allow the function to determine it by fitting the data.  In cases where the Ord plot heuristic recommends the Poissonness plot, it also provides this parameter, which I provide to the distplot function for the lower right plot.  Thus, while both the lower right and lower left plots are Poissonness plots in this case, they are generally based on different distribution parameters.  In the particular example shown here – constructed from the “number of times pregnant” variable in the Pima Indians diabetes dataset that I have discussed in several previous posts (available from the UCI Machine Learning Repository) – the Ord plot heuristic recommends the negative binomial distribution.  Comparing the Poissonness and negative-binomialness plots in the bottom row of the above plot array, it does appear that the negative binomial distribution fits the data better.

Finally, before examining the code for the CountSummary procedure, it is worth noting that the vcd package’s implementation of the Ord_plot and Ord_estimate procedures can generate four different distributional recommendations: the Poisson and negative-binomial distributions discussed here, along with the binomial distribution and the much less well-known log-series distribution.  The distplot procedure is flexible enough to generate plots for the first three of these distributions, but not the fourth, so in cases where the Ord plot heuristic recommends this last distribution, the CountSummary procedure displays the recommended distribution and parameter, but displays a warning message that no distribution plot is available for this case in the lower right plot position.

The code for the CountSummary procedure looks like this:

CountSummary <- function(xCount,TitleString){

  #

  #  Initial setup

  #

  library(vcd)

  grid.newpage()

  #

  #  Set up 2x2 array of plots

  #

  pushViewport(viewport(layout=grid.layout(nrow=2,ncol=2)))

  #

  #  Generate the plots:

  #

  #  1 - upper left = Ord plot

  #

  pushViewport(viewport(layout.pos.row=1,layout.pos.col=1))

  OrdLine = Ord_plot(xCount, newpage = FALSE, pop=FALSE, legend=FALSE)

  OrdType = Ord_estimate(OrdLine)

  popViewport()

  #

  #  2 - upper right = text summary

  #

  OrdTypeText = paste("Type = ",OrdType$type,sep=" ")

  if (OrdType$type == "poisson"){

    OrdPar = "Lambda = "

  }

  else if ((OrdType$type == "nbinomial")|(OrdType$type == "nbinomial")){

    OrdPar = "Prob = "

  } 

  else if (OrdType$type == "log-series"){

    OrdPar = "Theta = "

  }

  else{

    OrdPar = "Parameter = "

  }

  OrdEstText = paste(OrdPar,round(OrdType$estimate,digits=3), sep=" ")

  TextSummary = paste("Ord plot heuristic results:",OrdTypeText,OrdEstText,sep="\n")

  pushViewport(viewport(layout.pos.row=1,layout.pos.col=2))

  grid.text(TitleString,y=2/3,gp=gpar(fontface="bold"))

  grid.text(TextSummary,y=1/3)

  popViewport()

  #

  #  3 - lower left = standard Poissonness plot

  #

  pushViewport(viewport(layout.pos.row=2,layout.pos.col=1))

  distplot(xCount, type="poisson",newpage=FALSE, pop=FALSE, legend = FALSE)

  popViewport() 

  #

  #  4 - lower right = plot suggested by Ord results

  #

  pushViewport(viewport(layout.pos.row=2,layout.pos.col=2))

  if (OrdType$type == "poisson"){

    distplot(xCount, type="poisson",lambda=OrdType$estimate, newpage=FALSE, pop=FALSE, legend=FALSE)

  }

  else if (OrdType$type == "nbinomial"){

    prob = OrdType$estimate

    size = 1/prob - 1

    distplot(xCount, type="nbinomial",size=size,newpage=FALSE, pop=FALSE, legend=FALSE)

  }

  else if (OrdType$type == "binomial"){

    distplot(xCount, type="binomial", newpage=FALSE, pop=FALSE, legend=FALSE)

  }

  else{

    Message = paste("No distribution plot","available","for this case",sep="\n")

    grid.text(Message)

  }

  popViewport()

  #

}

This procedure is a function called with two arguments: the sequence of count values, xCounts, and TitleString, a text string that is displayed in the upper right text box in the plot array, along with the recommendations from the Ord plot heuristic.  When called, the function first loads the vcd library to make the Ord_plot, Ord_estimate, and distplot functions available for use, and it executes the grid.newpage() command to clear the display.  (Note that we don’t have to include “library(grid)” here to load the grid package, since loading the vcd package automatically does this.)   As in the previous example, the first “pushViewport” command creates the two-by-two plot array, and this is again followed by four code segments, one to generate each of the four displays in this array.  The first of these segments invokes the Ord_plot and Ord_estimate commands as discussed above, first to generate the upper left plot (a side-effect of the Ord_plot command) and second, to obtain the Ord plot heuristic recommendations, to be used in structuring the rest of the display.  The second segment creates a text display as in the first example considered here, but the structure of this display depends on the Ord plot heuristic results (i.e., the names of the parameters for the four possible recommended distributions are different, and the logic in this code block matches the display text to this distribution).  As noted in the preceding discussion, the third plot (lower left) is the Poissonness plot generated by the distplot function from the vcd package.  In this case, the function is called only specifying ‘type = “poisson”’ without the optional distribution parameter lambda, which is obtained by fitting the data.  The final element of this plot array, in the lower right, is also generated via a call to the distplot function, but here, the results from the Ord plot heuristic are used to specify both the type parameter and any optional or required shape parameters for the distribution.  As with the displayed text, simple if-then-else logic is required here to match the plot generated with the Ord plot heuristic recommendations.

Finally, it is important to note that in all of the calls made to Ord_plot or distplot in the CountSummary procedure, the parameters newpage, pop, and legend, are all specified as FALSE.  Specifying “newpage = FALSE” prevents these vcd plot commands from clearing the display page and erasing everything we have done so far.  Similarly, specifying “pop = FALSE” allows us to continue working in the current plot window until we notify the grid graphics system that we are done with it by issuing our own “popViewport()” command.  Specifying “legend = FALSE” tells Ord_plot and distplot not to write the default informational legend on each plot.  This is important here because, given the relatively small size of the plots generated in this two-by-two array, including the default legends would obscure important details.

Base versus grid graphics

In a comment in response to my latest post, Robert Young took issue with my characterization of grid as an R graphics package. Perhaps grid is better described as a “graphics support package,” but my primary point – and the main point of this post – is that to generate the display you want, it is sometimes necessary to use commands from this package. In my case, the necessity to learn something about grid graphics came as the result of my attempt to implement the CountSummary procedure to be included in the ExploringData package that I am developing. CountSummary is a graphical summary procedure for count data, based on Poissonness plots, negative binomialness plots, and Ord plots, all discussed in Chapter 8 of Exploring Data in Engineering, the Sciences and Medicine. My original idea was to implement these plots myself, but then I discovered that all three were already available in the vcd package. One of the great things about R is that you are encouraged to build on what already exists, so using the vcd implementations seemed like a no-brainer. Unfortunately, my first attempt at creating a two-by-two array of plots from the vcd package failed, and I didn’t understand why. The reason turned out to be that I was attempting to mix the base graphics command “par(mfrow=c(2,2))” that sets up a two-by-two array with varous plotting commands from vcd, which are based on grid graphics. Because these two graphics systems don’t play well together, I didn’t get the results I wanted. In the end, however, by learning a little about the grid package and its commands, I was able to generate my two-by-two plot array without a great deal of difficulty. Since grid graphics isn’t even mentioned in my favorite R reference book (Michael Crawley’s The R Book), I wanted to say a little here about what the grid package is and why you might need to know something about it. To do this, I will describe the ideas that went into the development of the CountSummary procedure and conclude this post with an example that shows what the output looks like. Next time, I will give a detailed discussion of the R code that generated these results. (For those wanting a preliminary view of what the code looks like, load the vcd package with the library command and run “examples(Ord_plot)” – in addition to generating the plots, this example displays the grid commands needed to construct the two-by-two array.)




Count variables – non-negative integer variables like the “number of times pregnant” (NPG) variable from the Pima Indians database described below – are often assumed to obey a Poisson distribution, in much the same way that continuous-valued variables are often assumed to obey a Gaussian (normal) distribution. Like this normality assumption for continuous variables, the Poisson assumption for count data is sometimes reasonable, but sometimes it isn’t. Normal quantile-quantile plots like those generated by the qqnorm command in base R or the qqPlot command from the car package are useful in informally assessing the reasonableness of the normality assumption for continuous data. Similarly, Poissonness plots are the corresponding graphical tool for informally evaluating the Poisson hypothesis for count data. The construction and interpretation of these plots is discussed in some detail in Chapters 8 and 9 of Exploring Data, but briefly, this plot constructs a variable called the Poissonness count metameter from the number of times each possible count value occurs in the data; if the data sequence conforms to the Poisson distribution, the points on this plot should fall approximately on a straight line. A simple R function that constructs Poissonness plots is available on the OUP companion website for the book, but an implementation that is both more conveniently available and more flexible is the distplot function in the vcd package, which also generates the negative binomialness plot discussed below.

The figure above is the Poissonness plot constructed using the distplot procedure from the vcd package for the NPG variable from the Pima Indians diabetes dataset mentioned above. I have discussed this dataset in previous posts and have used it as the basis for several examples in Exploring Data. It is available from the UCI Machine Learning Repository and it has been incorporated in various forms as an example dataset in a number of R packages, including a cleaned-up version in the MASS package (dataset Pima.tr). The full version considered here contains nine characteristics for 768 female members of the Pima Indian tribe, including their age, medical characteristics like diastolic blood pressure, and the number of times each woman has been pregnant. If this NPG count sequence obeyed the Poisson distribution, the points in the above plot would fall approximately on the reference line included there. The fact that these points do not conform well to this line – note, in particular, the departure at the lower left end of the plot where most of the counts occur – calls the Poisson working assumption into question.



A fundamental feature of the Poisson distribution is that it is defined by a single parameter that determines all distributional characteristics, including both the mean and the variance. In fact, a key characteristic of the Poisson distribution is that the variance is equal to the mean. This constraint is not satisfied by all count data sequences we encounter, however, and these deviations are important enough to receive special designations: integer sequences whose variance is larger than their mean are commonly called overdispersed, while those whose variance is smaller than their mean are commonly called underdispersed. In practice, overdispersion seems to occur more frequently, and a popular distributional alternative for overdispersed sequences is the negative binomial distribution. This distribution is defined by two parameters and it is capable of matching both the mean and variance of arbitrary overdispersed count data sequences. For a detailed discussion of this distribution, refer to Chapter 3 of Exploring Data.



Like the Poisson distribution, it is possible to evaluate the reasonableness of the negative binomial distribution graphically, via the negative binomialness plot. Like the Poissonness plot, this plot is based on a quantity called the negative binomialness metameter, computed from the number of times each count value occurs, plotted against those count values. To construct this plot, it is necessary to specify a numerical value for the distribution’s second parameter (the size parameter in the distplot command, corresponding to the r parameter in the discussion of this distribution given in Chapter 8 of Exploring Data). This can be done in several different ways, including the specification of trial values, the approach taken in the negative binomialness plot procedure that is available from the OUP companion website. This option is also available with the distplot command from the vcd package: to obtain a negative binomialness plot, specify the type parameter as “nbinomial” and, if a fixed size parameter is desired, it is specified by giving a numerical value for the size parameter in the distplot function call. Alternatively, if this parameter is not specified, the distplot procedure will estimate it via the method of maximum likelihood, an extremely useful feature, although it is important to note that this estimation process can be time-consuming, especially for long data sequences. Finally, a third approach that can be adopted is to use the Ord plot described next to obtain an estimate of this parameter based on a simple heuristic. In addition, this heuristic suggests which of these two candidate distributions – the Poisson or the negative binomial – is more appropriate for the data sequence.



Like the Poissonness plot, the Ord plot computes a simple derived quantity from the original count data sequence – specifically, the frequency ratio, defined for each count value as that value multiplied by the ratio of the number of times it occurs to the number of times the next smaller count occurs – and plots this versus the counts. If the data sequence obeys the negative binomial distribution, these points should conform reasonably well to a line with positive slope, and this slope can be used to determine the size parameter for the distribution. Conversely, if the Poisson distribution is appropriate, the best fit reference line for the Ord plot should have zero slope. In addition, Ord plots can also be used to suggest two additional discrete distributions (specifically, the binomial distribution and the log-series distribution), and the vcd package provides dataset examples to illustrate all four of these cases.



For my CountSummary procedure, I decided to construct a two-by-two array with the following four components. First, in the upper left, I used the Ord_plot command in vcd to generate an Ord plot. This command returns the intercept and slope parameters for the reference line in the plot, and the Ord_estimate command can then be used to convert these values into a type specification and an estimate of the distribution parameter needed to construct the appropriate discrete distribution plot. I will discuss these results in more detail in my next post, but for the case of the NPG count sequence considered here, the Ord plot results suggest the negative binomial distribution as the most appropriate choice, returning a parameter prob, from which the size parameter required to generate the negative binomialness plot may be generated (specifically, size = 1/prob – 1). The upper right quadrant of this display gives a text summary identifying the variable being characterized and listing the Ord plot recommendations and parameter estimate. Since the Poisson distribution is “the default” assumption for count data, the lower left plot shows a Poissonness plot for the data sequence, while the lower right plot is the “distribution-ness plot” for the distribution recommended by the Ord plot results. The results obtained by the CountSummary procedure for the NPG sequence are shown below. Next time, I will present the code used to generate this plot.

Graphical insights from the 2012 UseR! Meeting

About this time last month, I attended the 2012 UseR! Meeting.  Now an annual event, this series of conferences started in Europe in 2004 as an every-other-year gathering that now seems to alternate between the U.S. and Europe.  This year’s meeting was held on the Vanderbilt University campus in Nashville, TN, and it was attended by about 500 R aficionados, ranging from beginners who have just learned about R to members of the original group of developers and the R Core Team that continues to maintain it.  Many different topics were discussed, but one given particular emphasis was data visualization, which forms the primary focus of this post.  For a more complete view of the range of topics discussed and who discussed them, the conference program is available as a PDF file that includes short abstracts of the talks.

All attendees were invited to present a Lightning Talk, and about 20 of us did.  The format is essentially the technical equivalent of the 50-yard dash: before the talk, you provide the organizers exactly 15 slides, each of which is displayed for 20 seconds.  The speaker’s challenge is first, to try to keep up with the slides, and second, to try to convey some useful information about each one.  For my Lightning Talk, I described the ExploringData R package that I am in the process of developing, as a companion to both this blog and my book, Exploring Data in Engineering, the Sciences, and Medicine.  The intent of the package is first, to make the R procedures and datasets from the OUP companion site for the book more readily accessible, and second, to provide some additional useful tools for exploratory data analysis, incorporating some of the extensions I have discussed in previous blog posts.

Originally, I had hoped to have the package complete by the time I gave my Lightning Talk, but in retrospect, it is just as well that the package is still in the development stage, because I picked up some extremely useful tips on what constitutes a good package at the meeting.  As a specific example, Hadley Wickham, Professor of Statistics at Rice University and the developer of the ggplot2 package (more on this later), gave a standing-room-only talk on package development, featuring the devtools package, something he developed to make the R package development process easier.  In addition, the CRC vendor display at the meeting gave me the opportunity to browse and purchase Paul Murrell’s book, R Graphics, which provides an extremely useful, detailed, and well-written treatment of the four different approaches to graphics in R that I will say a bit more about below.

Because I am still deciding what to include in the ExploringData package, one of the most valuable sessions for me was the invited talk by Di Cook, Professor of Statistics at Iowa State University, who emphasized the importance of meaningful graphical displays in understanding the contents of a dataset, particularly if it is new to us.  One of her key points – illustrated with examples from some extremely standard R packages – was that the “examples” associated with datasets included in R packages often fail to include any such graphical visualization, and even for those that do, the displays are often too cryptic to be informative.  While this point is obvious enough in retrospect, it is one that I – along with a lot of other people, evidently – had not thought about previously.  As a consequence, I am now giving careful thought to the design of informative display examples for each of the datasets I will include in the ExploringData package.

As I mentioned above, there are (at least) four fundamental approaches to doing graphics in R.  The one that most of us first encounter – the one we use by default every time we issue a “plot” command – is called base graphics, and it is included in base R to support a wide range of useful data visualization procedures, including scatter plots, boxplots, histograms, and a variety of other common displays.  The other three approaches to graphics – grid graphics, lattice graphics, and ggplot2 – all offer more advanced features than what is typically available in base graphics, but they are, most unfortunately, incompatible in a number of ways with base graphics.  I discovered this the hard way when I was preparing one of the procedures for the ExploringData package (the CountSummary procedure, which I will describe and demonstrate in my next post).  Specifically, the vcd package includes implementations of Poissonness plots, negative binomialness plots, and Ord plots, all discussed in Exploring Data, and I wanted to take advantage of these implementations in building a simple graphical summary display for count data.  In base graphics, to generate a two-by-two array of plots, you simply specify “par(mfrow=c(2,2))” and then generate each individual plot using standard plot commands.  When I tried this with the plots generated by the vcd package, I didn’t get what I wanted – for the most part, it appeared that the “par(mfrow=c(2,2))” command was simply being ignored, and when it wasn’t, multiple plots were piled up on top of each other.  It turns out that the vcd package uses grid graphics, which has a fundamentally different syntax: it’s more complicated, but in the end, it does provide a wider range of display options.  Ultimately, I was able to generate the display I wanted, although this required some digging, since grid graphics aren’t really discussed much in my standard R reference books.  For example, The R Book by Michael J. Crawley covers an extremely wide range of useful topics, but the only mentions of “grid” in the index refer to the generation of grid lines (e.g., the base graphics command “grid” generates grid lines on a base R plot, which is not based on grid graphics). 

Often, grid graphics are mentioned in passing in introductory descriptions of trellis (lattice) graphics, since the lattice package is based on grid graphics.  This package is discussed in The R Book, and I have used it occasionally because it does support things like violin plots that are not part of base graphics.  To date, I haven’t used it much because I find the syntax much more complicated, but I plan to look further into it, since it does appear to have a lot more capability than base graphics do.  Also, Murrell’s R Graphics book devotes a chapter to trellis graphics and the lattice package, which goes well beyond the treatments given in my other R references, and this provides me further motivation to learn more.  The fourth approach to R graphics – Hadley Wickham’s ggplot2 package – was much discussed at the UseR! Meeting, appearing both in examples presented in various authors’ talks and as components for more complex and specialized graphics packages.  I have not yet used ggplot2, but I intend to try it out, since it appears from some of the examples that this package can generate an extremely wide range of data visualizations, with simple types comparable to what is found in base graphics often available as defaults.  Like the lattice package, ggplot2 is also based on grid graphics, making it, too, incompatible with base graphics.  Again, the fact that Murrell’s book devotes a chapter to this package should also be quite helpful in learning when and how to make the best use of it.

This year’s UseR! Meeting was the second one I have attended – I also went to the 2010 meeting in Gaithersburg, MD, held at the National Institute of Standards and Technology (NIST).  Both have been fabulous meetings, and I fully expect future meetings to be as good: next year’s UseR! meeting is scheduled to be held in Spain and I’m not sure I will be able to attend, but I would love to.  In any case, if you can get there, I highly recommend it, based on my experiences so far.

Classifying the UCI mushrooms

In my last post, I considered the shifts in two interestingness measures as possible tools for selecting variables in classification problems.  Specifically, I considered the Gini and Shannon interestingness measures applied to the 22 categorical mushroom characteristics from the UCI mushroom dataset.  The proposed variable selection strategy was to compare these values when computed from only edible mushrooms or only poisonous mushrooms.  The rationale was that variables whose interestingness measures changed a lot between these two subsets might be predictive of mushroom edibility.  In this post, I examine this question a little more systematically, primarily to illustrate the mechanics of setting up classification problems and evaluating their results.

More specifically, the classification problem I consider here is that of building and comparing models that predicts mushroom edibility, each one based on a different mushroom characteristic.  In practice, you would generally consider more than one characteristic as the basis for prediction, but here, I want to use standard classification tools to provide a basis for comparing the predictabilities of each of the potentially promising mushroom characteristics identified in my last post.  In doing this, I also want to highlight three aspects of classification problems: first, the utility of randomly splitting the available data into subsets before undertaking the analysis, second, the fact that we have many different options in building classifiers, and third, one approach to assessing classification results.

One of the extremely useful ideas emphasized in the machine learning literature is the utility of randomly partitioning our dataset into three parts: one used to fit whatever prediction model we are interested in building, another used to perform intermediate fit comparisons (e.g., compare the performance of models based on different predictor variables), and a third that is saved for a final performance assessment.  The reasoning behind this partitioning is that if we allow our prediction model to become too complex, we run the risk of overfitting, or predicting some of the random details in our dataset, resulting in a model that does not perform well on other, similar datasets.  This is an important practical problem that I illustrate with an extreme example in Chapter 1 of Exploring Data in Engineering, the Sciences, and Medicine.  There, a sequence of seven monotonically-decaying observations is fit to a sixth-degree polynomial that exactly predicts the original seven observations, but which exhibits horrible interpolation and extrapolation behavior.  The point here is that we need a practical means of protecting ourselves against building models that are too specific to the dataset at hand, and the partitioning strategy just described provides a simple way of doing this.  That is, once we partition the data, we can fit our prediction model to the first subset and then evaluate its performance with respect to the second subset: because these subsets were generated by randomly sampling the original dataset, their general character is the same, so a “good” prediction model built from the first subset should give “reasonable” predictions for the second subset.  The reason for saving out a third data subset – not used at all until the final evaluation of our model – is that model-building is typically an iterative procedure, so we are likely to cycle repeatedly between the first and second subsets.  For the final model evaluation, it is desirable to have a dataset available that hasn’t been used at all.

Generating this three-way split in R is fairly easy.  As with many tasks, this can be done in more than one way, but the following procedure is fairly straightforward and only makes use of procedures available in base R:

RandomThreeWay.proc <- function(df, probs = c(35,35,30), iseed = 101){

  #

  set.seed(iseed)

  n = nrow(df)

  #

  u = runif(n)

  #

  nprobs = probs/sum(probs)

  brks = c(0,cumsum(nprobs))

  Subgroup = cut(u, breaks=brks, labels=c("A","B","C"), include.lowest=TRUE)

  #

  Subgroup

  #

}

This function is called with three parameters: the data frame that we wish to partition for our analysis, a vector of the relative sizes of our three partitions, and a seed for the random number generator.  In the implementation shown here, the vector of relative sizes is given the default values 35%/35%/30%, but any relative size partitioning can be specified.  The result returned by this procedure is the character vector Subgroup, which has the values “A”, “B”, or “C”, corresponding to the three desired partitions of the dataset.  The first line of this procedure sets the seed for the uniform random number generator used in the third line, and the second line specifies how many random numbers to generate (i.e., one for each data record in the data frame).  The basic idea here is to generate uniform random numbers on the interval [0,1] and then assign subgroups depending on whether this value falls into the interval between 0 and 0.35, 0.35 to 0.70, or 0.70 to 1.00.  The runif function generates the required random numbers, the cumsum function is used to generate the cumulative breakpoints from the normalized probabilities, and the cut function is used to group the uniform random numbers using these break points.

In the specific example considered here, I use logistic regression as my classifier, although many, many other classification procedures are available in R, including a wide range of decision tree-based models, random forest models, boosted tree models, naïve Bayes classifiers, and support vector machines, to name only a few.  (For a more complete list, refer to the CRAN task view on Machine Learning and Statistical Learning).  Here, I construct and compare six logistic regression models, each constructed to predict the probability that a mushroom is poisonous from one of the six mushroom characteristics identified in my previous post: GillSize, StalkShape, CapSurf, Bruises, GillSpace, and Pop.  In each case, I extract the records for subset “A” of the UCI mushroom dataset, as described above, and use the base R procedure glm to construct a logistic regression model.  Because the model evaluation procedure (somers2, described below) that I use here requires a binary response coded as 0 or 1, it is simplest to construct a data frame with this response explicitly, along with the prediction covariate of interest.  The following code does this for the first predictor (GillSize):

EorP = UCImushroom.frame$EorP

PoisonBinary = rep(0,length(EorP))

PoisonIndx = which(EorP = = "p")

PoisonBinary[PoisonIndx] = 1

FirstFrame = data.frame(PoisonBinary = PoisonBinary, Covar = UCImushroom.frame$GillSize)

In particular, this code constructs a two-column data frame that contains the binary response variable PoisonBinary that is equal to 1 whenever EorP is “p” and 0 whenever this variable is “e”, and the prediction covariate Covar, which is here “GillSize”.  Given this data frame, I then apply the following code to randomly partition this data frame into subsets A, B, and C, and I invoke the built-in glm procedure to fit a logistic regression model:

Subset = RandomThreeWay.proc(FirstFrame)

IndxA = which(Subset = = "A")

LogisticModel = glm(PoisonBinary ~ Covar, data = FirstFrame, subset = IndxA, family=binomial())

Note that here I have specified the model form using the R formula construction “PoisonBinary ~ Covar”, I have used the subset argument of the glm procedure to specify that I only want to fit the model to subset A, and I have specified “family = binomial()” to request a logistic regression model.  Once I have this model, I evaluate it using the concordance index C available from the somers2 function in the R package Hmisc.  This value corresponds to the area under the ROC curve and is a measure of agreement between the predictions of the logistic regression model and the actual binary response.  As discussed above, I want to do this evaluation for subset B to avoid an over-optimistic view of the model’s performance due to overfitting of subset A.  To do this, I need the model predictions from subset B, which I obtain with the built-in predict procedure:

IndxB = which(Subset = = "B")

PredPoisonProb = predict(LogisticModel, newdata = FirstFrame[IndxB,], type="response")

ObsPoisonBinary = FirstFrame$PoisonBinary[IndxB]

In addition, I have created the variable ObsPoissonBinary, the sequence of binary responses from subset B, which I will use in calling the somers2 function:

library(Hmisc)

somers2(PredPoisonProb, ObsPoisonBinary)

           C                     Dxy                    n                         Missing

   0.7375031    0.4750063 2858.0000000    0.0000000

The results shown here include the concordance index C, an alternative (and fully equivalent) measure called Somers’ D (from which the procedure gets its name), the number of records in the dataset (here, in subset B), and the number of missing records (here, none).  The concordance index C is a number that varies between 0 and 1, with values between 0.5 and 1.0 meaning that the predictions are better than random guessing, and values less than 0.5 indicating performance so poor that it is actually worse than random guessing.  Here, the value of approximately 0.738 suggests that GillSize is a reasonable predictor of mushroom edibility, at least for mushrooms like those characterized in the UCI mushroom dataset. 

Repeating this process for all six of the mushroom characteristics identified as potentially predictive by the interestingness change analysis I discussed last time leads to the following results:

           

            Pop:                 C = 0.753        (6 levels)

            Bruises:            C = 0.740        (2 levels)

            GillSize:            C = 0.738        (2 levels)

            GillSpace:         C = 0.635        (2 levels)

            CapSurf:           C = 0.595        (4 levels)

            StalkShape:      C = 0.550        (2 levels)

These results leave open the questions of whether other mushroom characteristics, not identified on the basis of their interestingness shifts, are in fact more predictive of edibility, or how much better the predictions can be if we use more than one prediction variable.  I will examine those questions in subsequent posts, using the ideas outlined here.  For now, it is enough to note that one advantage of the approach described here, relative to that using odds ratios for selected covariates discussed last time, is that this approach can be used to assess the potential prediction power of categorical variables with arbitrary numbers of levels, while the odds ratio approach is limited to two-level predictors.

Interestingness comparisons

In three previous posts (April 3, 2011,  April 12, 2011,and May 21, 2011), I have discussed interestingness measures, which characterize the distributional heterogeneity of categorical variables.  Four specific measures are discussed in Chapter 3 of Exploring Data in Engineering, the Sciences and Medicine: the Bray measure, the Gini measure, the Shannon measure, and the Simpson measure.  All four of these measures vary from 0 to 1 in value, exhibiting their minimum value when all levels of the variable are equally represented, and exhibiting their maximum value when the variable is completely concentrated on a single one of its several possible levels.  Intermediate values correspond to variables that are more or less homogeneously distributed: more homogeneous for smaller values of the measure, and less homogeneous for larger values.  One of the points I noted in my first post on this topic was that the different measures exhibit different behavior for the intermediate cases, reflecting different inherent sensitivities to the various ways in which a variable can be “more homogeneous” or “less homogeneous.”  This post examines changes in interestingness measures as a potential exploratory analysis tool for selecting categorical predictors of some binary response.  In fact, I examined the same question from a different perspective in my April 12 post noted above: the primary difference is that there, the characterization I considered generates a single graph for each variable, with the number of points on the graph corresponding to the number of levels of the variable.  Here, I examine a characterization that represents each variable as a single point on the graph, allowing us to consider all variables simultaneously.

As a reminder of how these measures behave, the figure above shows a plot of the normalized Gini measure versus the normalized Shannon measure for the 23 categorical variables included in the mushroom dataset from the UCI Machine Learning Repository.  As I have noted in several previous posts that have discussed  this dataset, it gives observable characteristics for 8,124 mushrooms and classifies each one as either edible or poisonous (the binary variable EorP).  The above plot illustrates the systematic difference between the normalized Shannon and Gini interestingness measures: there, each point represents one of the 23 variables in the dataset, with the horizontal axis representing the Shannon measure computed for the variable and the vertical axis rperesenting the corresponding Gini measure.  The plot shows that the Gini measure is consistently larger than the Shannon measure, since all points lie above the equality reference line in this plot except for the single point at the origin.  This point corresponds to the variable VeilType, which only exhibits a single value in this dataset, meaning that both the Gini and Shannon measures are inherently ill-defined; consequently, they are given the default value of zero here, consistent with the general interpretation of these measures: if a variable only assumes a single value, it seems reasonable to consider it “completely homogeneous.”

Because edible and poisonous mushrooms are fairly evenly represented in this dataset (51.8% edible versus 48.2% poisonous), it has been widely used as one of several benchmarks for evaluating classification algorithms.  In particular, given the other mushroom characteristics, the fundamental classification question is how well can we predict whether each mushroom is poisonous or edible.  In this post and a subsequent follow-up post, I consider a closely related question: can differences in a variable’s interestingness measure between the edible subset and the poisonous subset be used to help us select prediction covariates for these classification algorithms?  In this post, I present some preliminary evidence to suggest that this may be the case, while in a subsequent post, I will put the question to the test by seeing how well the covariates suggested by this analysis actually predict edibility.

The specific idea I examine here is the following: given an interestingness measure and a mushroom characteristic, compute this measure for the chosen characteristic, applied the edible and poisonous mushrooms separately.  If these numbers are very different, this suggests that the distribution of levels is different for edible and poisonous mushrooms, further suggesting that this variable may be a useful predictor of edibility.  To turn this idea into a data analysis tool, it is necessary to define what we mean by “very different,” and this can be done in more than one way.  Here, I consider two possibilities.  The first is what I call the “normalized difference,” defined as the difference of the two interestingness measures divided by their sum.  Since both interestingness measures lie between 0 and 1, it is not difficult to show that this normalized difference lies between -1 and +1.  As a specific application of this idea, consider the plot below, which shows the normalized difference in the Gini measure between the poisonous mushrooms and the edible mushrooms (the normalized Gini shift) plotted against the corresponding difference for the Shannon measure (the normalized Shannon shift).  In addition, this plot shows an equality reference line, and the fact that the points consistently lie between this line and the horizontal axis shows that the normalized Gini shift is consistently smaller in magnitude than the normalized Shannon shift.  This suggests that the normalized Shannon measure may be more sensitive to distributional differences between edible and poisonous mushrooms.

The next figure, below, shows a re-drawn version of the above plot, with the equality reference line removed and replaced by four other reference lines.  The vertical dashed lines correspond to the outlier detection limits obtained by the Hampel identifier with threshold value t = 2 (see Chapter 7 of Exploring Data for a detailed discussion of this procedure), computed from the normalized Shannon shift values, while the horizontal dashed lines represent the corresponding limits computed from the normalized Gini shift values.  Points falling outside these limits represent variables whose changes in both Gini measure and Shannon measure are “unusually large” according to the Hampel identifier criteria used here.  These points are represented as solid circles, while those not detected as “unusual” by the Hampel identifier are represented as open circles.  The idea proposed here – to be investigated in a future post – is that these outlying variables may be useful in predicting mushroom edibility.

More specifically, the five solid circles in the above plot correspond to the following mushroom characteristics.  The two points in the lower left corner of the plot – exhibiting almost the most negative normalized Shannon shift possible – correspond to GillSize and StalkShape, two binary variables.  As I discussed in a previous post (May 7, 2011) and I discuss further in Chapter 13 of Exploring Data, an extremely useful measure of association between two binary variables (e.g., between GillSize and edibility) is the odds ratio.  An examination of the odds ratios for these two variables suggest that both should be at least somewhat predictive of edibility: the odds ratio between GillSize and edibility is 0.056, suggesting a very strong association (specifically, a GillSize value of “n” for “narrow” is most commonly associated with poisonous mushrooms in the UCI mushroom dataset), while the odds ratio between StalkShape and edibility is less extreme at 1.511, but still different enough from the neutral value of 1 to be suggestive of a clear association between these variables (a StalkShape value of “t” is more strongly associated with edible mushrooms than the alternative value of “e”).  The solid circle in the upper right of this plot corresponds to the variable CapSurf, which has four levels and whose distributional homogeneity appears to change quite substantially, according to both the Gini and Shannon measures.  Because this variable has more than two levels, it is not possible to characterize its association in terms of its odds ratio relative to edibility.  Finally, the cluster of three points in the upper right, just barely above the upper horizontal dashed line, correspond to the binary variables Bruises and GillSpace, and the six-level variable Pop.  Both of these binary variables exhibit very large odds ratios with respect to edibility (9.97 and 13.55 for Bruises and GillSpace, respectively), again suggesting that these variables may be highly predictive of edibility. 

The prevalence of binary variables in these results is noteworthy, and it reflects the fact that distributional shifts for binary variables can only occur in one way (i.e., the relative frequency of either fixed level can either increase or decrease).  Thus, large shifts in either interestingness measure should correspond to significant odds ratios with respect to the binary response variable, and this is seen to be the case here.  The situation is more complicated when a variable exhibits more than two levels, since the distribution of these levels can change in many ways between the two binary response values.  An important advantage of techniques like the the interestingness shift analysis described here is that they are not restricted to binary characteristics, as odds ratio characterizations are.

The second approach I consider for measuring the shift in interestingness between edible and poisonous mushrooms is what I call the “marginal measure,” corresponding to the difference in either the Gini or the Shannon measure between poisonous and edible mushrooms, divided by the original measure for the complete dataset.  An important difference between the marginal measure and the normalized measure is that the marginal measure is not bounded to lie between -1 and +1, as is evident in the plot below.  This plot shows the marginal Gini shift against the marginal Shannon shift for the mushroom characteristics, in the same format as the plot above.  Here, only four points are flagged as outliers, corresponding to the four binary variables identified above from the normalized shift plot: Bruises (the point in the extreme upper right), GillSpace (the point just barely in the upper right quadrant), and GillSize and StalkShape (the two points in the extreme lower left).  However, if we lower the Hampel identifier threshold from t = 2 to t = 1.5, we again identify CapSurf and Pop as potentially influential variables.

This last observation suggests an alternative interpretation approach that may be worth exploring.  Specifically, both of the two previous plots give clear visual evidence of “cluster structure,” and the Hampel identifier does extract some or all of this structure from the plot, but only if we apply a sufficiently judicious tuning to the threshold parameter.  A possible alternative would be to apply cluster analysis procedures, and this will be the subject of one or more subsequent posts.  In particular, there are many different clustering algorithms that could be applied to this problem, and the results are likely to be quite different.  The key practical question is which ones – if any – lead to useful ways of grouping these mushroom characteristics.  Subsequent posts will examine this question further from several different perspectives.

David Olive’s median confidence interval

As I have discussed in a number of previous posts, the median represents a well-known and widely-used estimate of the “center” of a data sequence.  Relative to the better-known mean, the primary advantage of the median is its much reduced outlier sensitivity.  This post briefly describes a simple confidence interval for the median that is discussed in a paper by David Olive, available on-line via the following link:

http://www.math.siu.edu/olive/ppmedci.pdf

As Olive notes in his paper and I further demonstrate in this post, an advantage of his confidence interval for the median is that it provides a simple, numerical way of identifying situations where the data values deserve a careful, graphical look.  In particular, he advocates comparing the traditional confidence interval for the mean with his confidence interval for the median: if these intervals are markedly different, it is worth investigating to understand why.  This strategy may be viewed as a particular instance of Collin Mallows’ “compute and compare” advice, discussed at the end of Chapter 7 of Exploring Data in Engineering, the Sciences, and Medicine.  The key idea here is that under “standard” working assumptions – i.e., distributional symmetry and approximate normality – the mean and the median should be approximately the same: if they are not, it probably means these working assumptions have been violated, due to outliers in the data, pronounced distributional asymmetry, or other less common phenomena like strongly multimodal data distributions or coarse quantization.  In the increasingly common case where we have a lot of numerical variables to consider, it may be undesirable or infeasible to examine them all graphically: numerical comparisons like the one described here may be automated and used to point us to subsets of variables that we really need to look at further.  In addition to describing this confidence interval estimator and illustrating it for three examples, this post also provides the R code to compute it. 

As a first example, the plot above shows the makeup flow rate dataset discussed in Exploring Data and available as the makeup dataset  (makeup.csv) from the book's companion website.  This plot shows 2,589 successive observations of the measured flow rate of a solvent recycle stream in an industrial manufacturing process.  In normal operation, this flow rate is just under 400 – in fact, the median flow rate is 393.86 – but this data record also includes measurements during time intervals when the process is either being shut down, is not running, or is being started back up, and during these periods the measured flow rates decrease toward zero, are approximately equal to zero, and increase from zero back to approximately 400, respectively.  Because of the presence of these anomalous segments in the data, the mean value is much smaller than the median: specifically, the mean is 315.46, actually serving as a practical dividing line between the normal operation segments (i.e., those data points that lie above the mean) and the shutdown segments (i.e., those data points that lie below the mean).  The dashed lines in this plot at 309.49 and 321.44 correspond to the classical 95% confidence interval for the mean, computed as described below.  In contrast, the dotted lines at 391.83 and 394.88 correspond to Olive’s 95% confidence interval for the median, also described below.  Before proceeding to a more detailed discussion of how these lines were determined, the three primary points to note from this figure are, first, that the two confidence intervals are very different (e.g., they do not overlap at all), second, that the mean confidence intervals are much wider than those for the median in this case, and third, that the median confidence interval lies well within the range of the normal operating data, while the mean confidence interval does not.  It is also worth noting that, if we simply remove the shutdown episodes from this dataset, the mean of this edited dataset is 397.7, a value that lies slightly above the upper 95% confidence interval for the median, but only slightly so (this and other data cleaning strategies for this dataset are discussed in some detail in Chapter 7 of Exploring Data).

Both the classical confidence interval for the mean and David Olive’s confidence interval for the median are based on the fact that these estimators are asymptotically normal: for a sufficiently large data sample, both the estimated mean and the estimated median approach the correct limits for the underlying data distribution, with a standard deviation that decreases inversely with the square root of the sample size.  Using this description directly would lead to confidence intervals based on the quantiles of the Gaussian distribution, but for small to moderate-sized samples, more accurate confidence intervals are obtained by replacing these Gaussian quantiles with those for the Student’s t-distribution with the appropriate number of degrees of freedom.  More specifically, for the mean, the confidence interval at a given level p is of the form:

            CI = (Mean – c_p SE, Mean + c_p SE),

where c_p is the constant derived from the Gaussian or Student’s t-distribution, and SE is the standard error of the mean, equal to the usual standard deviation estimate divided by the square root of the number of data points.  (For a more detailed discussion of the math behind these results, refer to either Chapter 9 of Exploring Data or to David Olive’s paper, available through the link given above.)  For the median, Olive provides a simple estimator for the standard error, described further in the next paragraph.  First, however, it is worth saying a little about the difference between the Gaussian and Student’s t-distribution in these results.  Probably the most commonly used confidence intervals are the 95% intervals – these are the confidence intervals shown in the plot above for the makeup flow rate data – which represent the interval that should contain the true distribution mean with probability at least 95%.  In the Gaussian case, the constant c_p for the 95% confidence interval is approximately 1.96, while for the Student’s t-distribution, this number depends on the degrees of freedom parameter.  In the case of the mean, the degrees of freedom is one less than the sample size, while for the median confidence intervals described below, this number is typically much smaller.  The difference between these distributions is that the c_p parameter decreases from a very large value for few degrees of freedom – e.g., the 95% parameter value is 12.71 for a single degree of freedom – to the Gaussian value (e.g., 1.96 for the 95% case) in the limit of infinite degrees of freedom.  Thus, using Student’s t-distribution instead of the Gaussian distribution results in wider confidence intervals, wider by the ratio of the Student’s t value for c_p to the Gaussian value.  The plot below shows this ratio for the 95% parameter c_p as the degree of freedom parameter varies between 5 and 200, with the dashed line corresponding to the Gaussian limit when this ratio is equal to 1.

The general structure of Olive’s confidence interval for the median is exactly analogous to that for the mean given above:

            CI = (Median – c_p SE, Median + c_p SE)

The key result of Olive’s paper is a simple estimator for the standard error SE, based on order statistics (i.e., rank-ordered data values like the minimum, median, and maximum).  Instead of describing these results mathematically, I have included an R procedure that computes the median, Olive’s standard error, the corresponding confidence intervals, and the classical results for the mean (again, for the mathematical details, refer to Olive’s paper; for a more detailed discussion of order statistics, refer to Chapter 6 of Exploring Data).  Specifically, the following R procedure is called with a vector y of numerical data values, and the default level of the resulting confidence interval is 95%, although this level can be changed by specifying an alternative value of alpha (this is 1 minus the confidence level, so alpha is 0.05 for the 95% case, 0.01 for 99%, etc.).

DOliveCIproc <- function(y, alpha = 0.05){

  #

  #  This procedure implements David Olive's simple

  #  median confidence interval, along with the standard

  #  confidence interval for the mean, for comparison

  #

  #  First, compute the median

  #

  n = length(y)

  ysort = sort(y)

  nhalf = floor(n/2)

  if (2*nhalf < n){

    #  n odd

    med = ysort[nhalf + 1]

  }

  else{

    # n even

    med = (ysort[nhalf] + ysort[nhalf+1])/2

  }

  #

  #  Next, compute Olive’s standard error for the median

  #

  Ln = nhalf - ceiling(sqrt(n/4))

  Un = n - Ln

  SE = 0.5*(ysort[Un] - ysort[Ln+1])

  #

  #  Compute the confidence interval based on Student’s t-distribution

  #  The degrees of freedom parameter p is discussed in Olive’s paper

  #

  p = Un - Ln - 1

  t = qt(p = 1 - alpha/2, df = p)

  medLCI = med - t * SE

  medUCI = med + t * SE

  #

  #  Next, compute the mean and its classical confidence interval

  #

  mu = mean(y)

  SEmu = sd(y)/sqrt(n)

  tmu = qt(p = 1 - alpha/2, df = n-1)

  muLCI = mu - tmu * SEmu

  muUCI = mu + tmu * SEmu

  #

  #  Finally, return a data frame with all of the results computed here

  #

  OutFrame = data.frame(Median = med, LCI = medLCI, UCI = medUCI,

                        Mean = mu, MeanLCI = muLCI, MeanUCI = muUCI,

                        N = n, dof = p, tmedian = t, tmean = tmu,

                        SEmedian = SE, SEmean = SEmu)

  OutFrame

}

Briefly, this procedure performs the following computations.  The first portion of the code computes the median, defined as the middle element of the rank-ordered list of samples if the number of samples n is odd, and the average of the two middle samples if n is even.  Note that the even/odd character of n is determined by using the floor function in R: floor(n/2) is the largest integer that does not exceed n/2.  Thus, if n is odd, the floor function rounds n/2 down to its integer part, so the product 2 * floor(n/2) is less than n, while if n is even, floor(n/2) is exactly equal to n/2, so this product is equal to n.  In addition, both the floor function and its opposite function ceiling are needed to compute the value Ln used in computing Olive’s standard error for the median.  The c_p values correspond to the parameters t and tmu that appear in this function, computed from the built-in R function qt (which returns quantiles of the t-distribution).  Note that for the median, the degrees of freedom supplied to this function is p, which tends to be much smaller than the degrees of freedom value n-1 for the mean confidence interval computed in the latter part of this function.

As a specific illustration of the results generated by this procedure, applying it to the makeup flow rate data sequence yields:

> DOliveCIproc(makeupflow)

    Median          LCI             UCI           Mean        MeanLCI    MeanUCI     N   dof         tmedian

1 393.3586   391.8338   394.8834   315.4609   309.4857   321.4361   2589  52   2.006647

     tmean      SEmedian     SEmean

1 1.960881   0.75987     3.047188

>

These results were used to construct the confidence interval lines in the makeup flow rate plot shown above.  In addition, note that these results also illustrate the point noted in the preceding discussion about the degrees of freedom used in constructing the Student’s t-based confidence intervals.  For the mean, the degrees of freedom is N-1, which is 2588 for this example, meaning that there is essentially no difference in this case between these confidence intervals and those based on the Gaussian limiting distribution.  In contrast, for the median, the degrees of freedom is only 52, giving a c_p value that is about 2.5% larger than the corresponding Gaussian case; for the next example, the degrees of freedom is only 16, making this parameter about 8% larger than the Gaussian limit.

One of the points I discussed in my last post was the instability of the median relative to the mean, a point I illustrated with the plot shown above.  This is a simulation-based dataset consisting of three parts: the first 100 points are narrowly distributed around the value +1, the 101^st point is exactly zero, and the last 100 points are narrowly distributed around the value -1.  As I noted last time, removing two points from either the first group or the last group can profoundly alter the median, while having very little effect on the mean.  The figure shown above includes, in addition to the data values, the 95% confidence intervals for both the mean (the dotted lines in the center of the plot) and the median (the heavy dashed lines at the top and bottom of the plot).  Here, the fact that the median confidence interval is enormously wider (by almost a factor of 13) than the mean confidence interval gives an indication of the instability of the median.  In fact, the data distribution in this example is strongly bimodal, corresponding to a case where order statistic-based estimators like the median and Olive’s standard error for it perform poorly, a point discussed in Chapter 7 of Exploring Data.

One of the other important cases where estimators based on order statistics can perform poorly is that of coarsely quantized data, such as temperatures recorded only to the nearest tenth of a degree.  The difficulty with these cases is that coarse quantization profoundly changes the nature of the data distribution.  Specifically, it is a standard result in statistics that the probability of any two samples drawn from a continuous distribution having exactly the same value is zero, but this is no longer true for discrete distributions (e.g., count data), and coarse quantization introduces an element of discreteness into the data distribution.  The above figure illustrates this point for a simple simulation-based example.  The upper left plot shows a random sample of size 200 drawn from a zero-mean, unit-variance Gaussian distribution, and the upper right plot shows the effects of quantizing this sample, rounding it to the nearest half-integer value.  The lower two plots are normal quantile-quantile plots generated by the R command qqPlot from the car package: in the lower left plot, almost all of the points fall within the 95% confidence interval around the normal reference line for this plot, while many of the points fall somewhat outside these confidence limits in the plot shown in the lower right.  The greatest difference, however, is in the “staircase” appearance of this lower right plot, reflecting the effects of the coarse quantization on this data sample: each “step” corresponds to a group of samples that have exactly the same value.

The influence of this quantization on Olive’s confidence interval for the median is profound: for the original Gaussian data sequence, the 95% confidence interval for the median is approximately (-0.222,0.124), compared with (-0.174,0.095) for the mean.  These results are consistent with our expectations: since the mean is the best possible location estimator for Gaussian data, it should give the narrower confidence interval, and it does.  For the quantized case, the 95% confidence interval for the mean is (-0.194, 0.079), fairly similar to that for the original data sequence, but the confidence interval for the median reduces to the single value zero.  This result represents an implosion of Olive’s standard error estimator for the median, exactly analogous to the behavior of the MADM scale estimate that I have discussed previously when a majority of the data values (i.e., more than 50% of them) are identical.  Here, the situation is more serious, since the MADM scale estimate does not implode for this example: the MADM scale for the original data sequence is 0.938, versus 0.741 for the quantized sequence.  The reason Olive’s standard error estimator is more prone to implosion in the face of coarse quantization is that it is based on a small subset of the original data sample.  In particular, the size of the subsample on which this estimator is based is p, the degrees of freedom for the t-distribution used in constructing the corresponding confidence interval, and this number is approximately the square root of the sample size.  Thus, for a sample of size 200 like the example considered here, MADM scale implosion requires just over half the sample to have the same value – 101 data points in this case – where Olive’s standard error estimator for the median can implode if 16 or more samples have the same value, and this is exactly what happens here: the median value is zero, and this value occurs 39 times in the quantized data sequence.

David Olive’s confidence interval for the median is easily computed and represents a useful adjunct to the median as a characterization of numerical variables.  As Olive advises, there is considerable advantage in computing and comparing both his median confidence interval and the corresponding standard confidence interval around the mean.  Although in the summary of his paper, Olive only mentions outliers as a potential cause of substantial differences between these two confidence intervals, this post has illustrated that disagreements can also arise from other causes, including light-tailed, bimodal, or coarsely quantized data, much like the situation with the MADM scale estimate versus the standard deviation.  In fact, as the last example discussed here illustrates, Olive’s standard error estimator for the median and the confidence intervals based on it can implode – exactly like the MADM scale estimate – in the face of coarsely quantized data.  In fact, the implosion problem for Olive’s median standard error estimator is potentially more severe, again as illustrated in the previous example.  Finally, it is worth noting that Olive’s paper also discusses confidence intervals for trimmed means.