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.

Gastwirth’s location estimator

The problem of outliers – data points that are substantially inconsistent with the majority of the other points in a dataset – arises frequently in the analysis of numerical data.  The practical importance of outliers lies in the fact that even a few of these points can badly distort the results of an otherwise reasonable data analysis.  This outlier-sensitivity problem is often particularly acute for classical data characterizations and analysis methods like means, standard deviations, and linear regression analysis.  As a consequence, a range of outlier-resistant methods have been developed for many different applications, and new methods continue to be developed.  For example, the R package robustbase that I have discussed in previous posts includes outlier-resistant methods for estimating location (i.e., outlier-resistant alternatives to the mean), estimating scale (outlier-resistant alternatives to the standard deviation), quantifying asymmetry (outlier-resistant alternatives to the skewness), and fitting regression models.  In Exploring Data in Engineering, the Sciences, and Medicine, I discuss a number of outlier-resistant methods for addressing some of these problems, including Gastwirth’s location estimator, an alternative to the mean that is the subject of this post.

The mean is the best-known location estimator, and it gives a useful assessment of the “typical” value of any numerical sequence that is reasonably symmetrically distributed and free of outliers.  The outlier-sensitivity of the mean is severe, however, which motivates the use of outlier-resistant alternatives like the median.  While the median is almost as well-known as the mean and extremely outlier-resistant, it can behave unexpectedly (i.e., “badly”) as a result of its non-smooth character.  This point is illustrated in Fig. 7.23 in Exploring Data, identical in character to the figure shown below (this figure is slightly different because it uses a different seed to generate the random numbers on which it is based).  Specifically, this plot shows a sequence of 201 data points, constructed as follows.  The first 100 points are normally distributed with mean 1 and standard deviation 0.1, the 101^st point is equal to zero, and points 102 through 201 are normally distributed with mean -1 and standard deviation 0.1.  Small changes in this dataset in the specific form of deleting points can result in very large changes in the computed median.  Specifically, in this example, the first 100 points lie between 0.768 and 1.185 and the last 100 points lie between -0.787 and -1.282; because the central data point lies between these two equal-sized groups, it defines the median, which is 0.  The mean is quite close to this value, at -0.004, but the situation changes dramatically if we omit either the first two or the last two points from this data sequence.  Specifically, the median value computed from points 1 through 199 is 0.768, while that computed from points 3 through 201 is -0.787.  In contrast, the mean values for these two modified sequences are 0.006 and -0.014.  Thus, although the median is much less sensitive than the mean to contamination from outliers, it is extremely sensitive to the 1% change made in this example for this particular dataset. 

The fact that the median is not “universally the best location estimator” provides a practical motivation for examining alternatives that are intermediate in behavior between the very smooth but very outlier-sensitive mean and the very outlier-insensitive but very non-smooth median.  Some of these alternatives were examined in detail in the book Robust Estimates of Location: Survey and Advances, by D.F. Andrews, P.J. Bickel, F.R. Hampel, P.J. Huber, W.H. Rogers, and J.W. Tukey, published by Princeton University Press in 1972 (according to the publisher's website, this book is out of print, but used copies are available through distributors like Amazon or Barnes and Noble).  The book summarizes the results of a year-long study of 68 different location estimators, including both the mean and the median.  The fundamental criteria for inclusion in this study were, first, that the estimators had to be computable from any given sequence of real numbers, and second, that they had to be both location and scale-invariant.  Specifically, if a given data sequence {x_k} yielded a result m, the scaled and shifted data sequence {Ax_k + b} should yield the result Am+b, for any numbers A and b.  The study was co-authored by six statistical researchers with differing opinions and points of view, but two of the authors – D.F. Andrews and F.R. Hampel – included the Gastwirth estimator (described in detail below) in their list of favorites.  For example, Hampel characterized this estimator as one of a small list of those that were “never bad at the distributions considered.”  Also, in contrast to many of the location estimators considered in the study, Gastwirth’s estimator does not require iterative computations, making it simpler to implement.

Specifically, Gastwirth’s location estimator is a weighted sum of three order statistics.  That is, to compute this estimator, we first sort the data sequence in ascending order.  Then, we take the values that are one-third of the way up this sequence (the 0.33 quantile), half way up the sequence (i.e., the median, or 0.50 quantile), and two-thirds of the way up the sequence (the 0.67 quantile).  Given these three values, we then form the weighted average, giving the central (median) value a weight of 40% and the two extreme values each a weight of 30%.  This is extremely easy to do in R, with the following code:

Gastwirth <- function(x,...){

  #

  ordstats = quantile(x, probs=c(1/3,1/2,2/3),...)

  wts = c(0.3,0.4,0.3)

  sum(wts*ordstats)

  #

}

The key part of this code is the first line, which computes the required order statistics (i.e., the quantiles 1/3, 1/2, and 2/3) using the built-in quantile function.  The first argument passed to this function is x, the vector of data values to be characterized, and the second argument (probs) defines the specific quantiles we wish to compute.  The ellipses in the Gastwirth procedure’s command line is passed to the quantile function; several parameters are possible (type “help(quantile)” in your R session for details), but one of the most useful is na.rm, a logical variable that specifies how missing data values are to be handled.  The default is “FALSE” and this causes the Gastwirth procedure to return the missing data value “NA” if any values of x are missing; the alternative “TRUE” computes the Gastwirth estimator from the non-missing values, giving a numerical result.  The three-element vector wts defines the quantile weights that define the Gastwirth estimator, which the final sum statement computes.

For the data example considered above, the Gastwirth estimator yields the location estimate -0.001 for the complete dataset, 0.308 for points 1 to 199 (vs. 0.768 for the median), and -0.317 for points 3 to 201 (vs. -0.787 for the median).  Thus, while it does not perform nearly as well as the mean for this example, it performs substantially better than the median. 

For the infinite-variance Cauchy distribution that I have discussed in several previous posts, the Gastwirth estimator performs similarly to the median, yielding a useful estimate of the center of the data distribution, in contrast to the mean, which doesn’t actually exist for this distribution (that is, the first moment does not exist for the Cauchy distribution).  Still, the distribution is symmetric about zero, so the median is well-defined, as is the Gastwirth estimator, and both should be zero for this distribution.  The above figure shows the results of applying these three estimators – the mean, the median, and Gastwirth’s estimator – to 1,000 independent random samples drawn from the Cauchy distribution.  Specifically, this figure gives a boxplot summary of these results, truncated to the range from -3 to 3 to show the range of variation of the median and Gastwirth estimator (without this restriction, the boxplot comparison would be fairly non-informative, since the mean values range from approximately -161 to 27,793, reflecting the fact that the mean is not a consistent location estimator for the Cauchy distribution).   To generate these results, the replicate function in R was used, followed by the apply function, as follows:

    RandomSampleFrame = replicate(1000, rt(n=200,df=1))
    BoxPlotVector = apply(RandomSampleFrame, MARGIN=2, Gastwirth)

The replicate function creates a data frame with the number of columns specified by the first argument (here, 1000), and each column generated by the R statement that appears as the second argument.  In this case, this second argument is the command rt, which generates a sequence of n statistically independent random numbers drawn from the Student’s t-distribution with the number of degrees of freedom specified by the df argument (here, this is 1, corresponding to the fact that the Cauchy distribution is the Student’s t-distribution with 1 degree of freedom).   Thus, RandomSampleFrame is a data frame with 200 rows and 1,000 columns, each of which may be regarded as a Cauchy-distributed random sample.  The apply function applies the function specified in the third argument (here, the Gastwirth procedure listed above) to the columns (MARGIN=2 specifies columns; MARGIN=1 would specify rows) of the data frame specified in the first argument.  The result is BoxPlotVector, a vector of 1,000 Gastwirth estimates, one for each random sample generated by the replicate function above.

At the other extreme, in the limit of infinite degrees of freedom, the Student’s t-distribution approaches a Gaussian limit.  The figure above shows the same comparison as before, except for the Gaussian distribution instead of the Cauchy distribution.  Here, the mean is the best possible location estimator and it clearly performs the best, but the point of this example is that Gastwirth’s location estimator performs better than the median.  In particular, the interquartile distance (i.e., the width of the “box” in each boxplot) for the mean is 0.094, it is 0.113 for the median, and it is 0.106 for Gastwirth’s estimator.

Another application area where very robust estimators like the median often perform poorly is that of bimodal distributions like the arc-sine distribution whose density is plotted above.  This distribution is a symmetric beta distribution, with both shape parameters equal to 0.5 (see Exploring Data, Sec. 4.5.1 for further discussion of this distribution).  Because it is symmetrically distributed on the interval from 0 to 1, the location parameter for this distribution is 0.5 and all three of the location estimators considered here yield values that are accurate on average, but with different levels of precision.  This point is shown in the figure below, which again provides boxplot comparisons for 1,000 random samples drawn from this distribution, each of length 200, for the mean, median, and Gastwirth location estimators.  As in the Gaussian case considered above, the mean performs best here, with an interquartile distance of 0.035, the median performs worst, with an interquartile distance of 0.077, and Gastwirth’s estimator is intermediate, with an interquartile distance of 0.060.

The point of this post has been to illustrate a location estimator with properties that are intermediate between those of the much better-known mean and median.  In particular, the results presented here for the Cauchy distribution show that Gastwirth’s estimator is intermediate in outlier sensitivity between the disastrously sensitive mean and the maximally insensitive median.  Similarly, the first example demonstrated that Gastwirth’s estimator is also intermediate in smoothness between the maximally smooth mean and the discontinuous median: the sensitivity of Gastwirth’s estimator to data editing in “swing-vote” examples like the one presented here is still undesirably large, but much better than that of the median.  Finally, the results presented here for the Gaussian and arc-sine distributions show that Gastwirth’s estimator is better-behaved for these distributions than the median.  Because it is extremely easy to implement in R, Gastwirth’s estimator seems worth knowing about.

Measuring associations between non-numeric variables

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

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

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

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

GKtau <- function(x,y){

  #

  #  First, compute the IxJ contingency table between x and y

  #

  Nij = table(x,y,useNA="ifany")

  #

  #  Next, convert this table into a joint probability estimate

  #

  PIij = Nij/sum(Nij)

  #

  #  Compute the marginal probability estimates

  #

  PIiPlus = apply(PIij,MARGIN=1,sum)

  PIPlusj = apply(PIij,MARGIN=2,sum)

  #

  #  Compute the marginal variation of y

  #

  Vy = 1 - sum(PIPlusj^2)

  #

  #  Compute the expected conditional variation of y given x

  #

  InnerSum = apply(PIij^2,MARGIN=1,sum)

  VyBarx = 1 - sum(InnerSum/PIiPlus)

  #

  #  Compute and return Goodman and Kruskal's tau measure

  #

  tau = (Vy - VyBarx)/Vy

  tau

}

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

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

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

            set.seed(123)

            u = runif(400)

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

            x = round(u,digits=1)

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

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

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

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

            c = letters[10*x+1]

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

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

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

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

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

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

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

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

Moving window filters and the pracma package

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

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

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

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

                        outlierMAD(x,k)

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

outlierMAD <- function (x, k){

    n <- length(x)

    y <- x

    ind <- c()

    L <- 1.4826

    t0 <- 3

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

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

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

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

            y[i] <- x0

            ind <- c(ind, i)

        }

    }

    list(y = y, ind = ind)

}

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

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

    n <- length(x)

    y <- x

    ind <- c()

    L <- 1.4826

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

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

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

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

            y[i] <- x0

            ind <- c(ind, i)

        }

    }

    list(y = y, ind = ind)

}

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

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

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

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

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

  #

  n = length(x)

  fl = 2*K+1

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

  if (dorder == 0){

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

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

  }

  else{

    y[1:K] = 0

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

  }

  y

}

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

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