Better Estimates from Binned Income Data: Interpolated CDFs and Mean-Matching

Researchers often estimate income statistics from summaries that report the number of incomes in bins such as $0 to 10,000, $10,001 to 20,000, …, $200,000+. Some analysts assign incomes to bin midpoints, but this treats income as discrete. Other analysts fit a continuous parametric distribution, but the distribution may not fit well. We fit nonparametric continuous distributions that reproduce the bin counts perfectly by interpolating the cumulative distribution function (CDF). We also show how both midpoints and interpolated CDFs can be constrained to reproduce the mean of income when it is known. We evaluate the methods in estimating the Gini coefficients of all 3,221 U.S. counties. Fitting parametric distributions is very slow. Fitting interpolated CDFs is much faster and slightly more accurate. Both interpolated CDFs and midpoints give dramatically better estimates if constrained to match a known mean. We have implemented interpolated CDFs in the “binsmooth” package for R. We have implemented the midpoint method in the “rpme” command for Stata. Both implementations can be constrained to match a known mean.


Introduction
Surveys often ask respondents to report income in brackets or bins, such as $0-10,000, $10,000-20,000,. . . ,$200,000+. Even in surveys where respondents report exact incomes, incomes may be binned before publication, either to protect privacy or to summarize the income distribution compactly with the number of incomes in each bin. Table 1 gives a binned summary of household incomes in Nantucket, the richest county in the US. Binning presents challenges to investigators who want to estimate simple summary statistics such as the mean, median, or standard deviation, or inequality statistics such as the Gini coefficient, the Theil index, the coefficient of variation, or the mean log deviation. Researchers have implemented several methods for calculating estimates from binned incomes.
The simplest and most popular approach is to assign each case to the midpoint of its bin-using a robust pseudo-midpoint for the top bin, whose upper bound is typically undefined (e.g., Table 1). The weakness of the midpoint approach is that it treats income as a discrete variable, but the method also has several strengths. The midpoint method is easy to implement and runs quickly. Midpoint estimates are also "bin-consistent"  in the sense that midpoint estimates will get arbitrarily close to their estimands if the bins are sufficiently numerous and narrow.
Another approach is to fit the bin counts to a continuous parametric distribution. Popular distributions include 2-, 3-, and 4-parameter distributions from the generalized beta family, which includes the Pareto, lognormal, Weibull, Dagum, and other distributions (McDonald and Ransom, 2008). One implementation fits up to 10 distributions and selects the one that fits best. An alternative is to use the AIC or BIC to calculated a weighted average of income statistics across several candidate distributions .
A strength of the parametric approach is that it treats income as continuous. A weakness is that even the best-fitting parametric distribution may not fit the bin counts particularly well. If the fit is poor, the parametric approach is not bin-consistent; that is, even with an infinite number of bins, each infinitesimally narrow, a parametric distribution may produce poor estimates if it is not a good fit to the underlying distribution of income.
A practical weakness of the parametric approach is that it is typically implemented using iterative methods which can be slow. The speed of a parametric fit may be acceptable if you fit a single distribution to a single binned dataset, but runtimes of hours are possible if you fit several distributions to thousands of binned datasets-such as every county or school district in the US. Other computational issues include nonconvergence and undefined estimates. These issues are rare but inevitable when you run thousands of binned datasets .
Neither approach-midpoint or parametric-is uniformly more accurate. With many bins, midpoint estimates are better because of bin-consistency, while with fewer than 8 bins, parametric estimates can be better because of their smoothness .
Empirically, the parametric and midpoint approaches produce similarly accurate estimates from typical US income data with 15 to 25 bins. Both methods typically estimate the Gini within a few percentage points of its true value. This is accurate enough for many purposes, but can lead to errors when estimating small differences or changes such as the 5% increase in the Gini of US family income that occurred between 1970and 1980.
A potential improvement is to fit binned incomes to a flexible nonparametric continuous density. Like parametric densities, a nonparametric density treats income as continuous.
Like the midpoint method, a nonparametric density can be bin-consistent and fit the bin counts as closely as we like.
Unfortunately, past nonparametric approaches have been disappointing. A nonparametric approach using kernel density estimation had substantial bias under some circumstances (Minoiu and Reddy, 2012). A nonparametric approach using a spline to model the log of the density (Kooperberg and Stone, 1992) had even greater bias (von Hippel et al., 2012), though it was not clear whether the bias came from the method or its software implementation.
In this paper, we implement and test a nonparametric continuous method that outperforms its predecessors in both speed and accuracy. The method, which we call CDF interpolation, simply connects points on the empirical cumulative distribution function (CDF). The method can connect the points using line segments or cubic splines. When cubic splines are used, the method is similar to "histospline" or "histopolation" methods which fit a spline to a histogram (Wahba, 1976;Morandi and Costantini, 1989). But histosplines are limited to histograms which have bins of equal width (Wang., 2015). CDF interpolation is more general approach which can handle income data where the bins have unequal width and the top bins has no upper bound (e.g., Table 1).
We have implemented CDF interpolation in our R package binsmooth (Hunter and Drown, 2016), which is available for download from the Comprehensive R Archive Network (CRAN).
Our results will show that statistics estimated with CDF interpolation are slightly more accurate than estimates obtained using midpoints or parametric distributions. In addition, CDF interpolation is much faster than parametric estimation, though not as fast as the midpoint method.
We also show that the differences between methods are dwarfed by the improvement we get if we constrain a method to match the grand mean of income, which the US Census Bureau often reports alongside the bin counts. If we constrain either the interpolated CDF or the midpoint method to match a known mean, we get dramatically better estimates of the Gini. Our binsmooth package can constrain an interpolated CDF to match a known mean, and our new version (2.0) of the rpme command for Stata can constrain the midpoint method to match a known mean as well.
In the rest of this paper, we define the midpoint, parametric, and interpolated CDF methods more precisely, then compare the accuracy of estimates in binned data summarizing household incomes within US counties. We also show how much estimates improve if we have the mean as well as the bin counts.

Binned data
A binned data set, such as Table 1, consists of counts n 1 , n 2 , . . . , n B specifying the number of cases in each of B bins. The total number of cases is T = n 1 + n 2 + · · · + n B . Each

A midpoint method
The oldest and simplest way to analyze binned data is the midpoint method, which within each bin b assigns incomes to the bin midpoint m b = (l b + u b )/2 (Heitjan1989 for a review).
Then statistics such as the Gini can be calculated by applying sample formulas to the midpoints m b weighted by the counts n b .
When the top bin has no upper bound, we must define a pseudo-midpoint for it. The traditional choice is µ B = l B α/(α − 1), which would define the arithmetic mean of the top bin if top-bin incomes followed a Pareto distribution with shape parameter α > 1 (Henson, 1967). The problem with this choice is that the arithmetic mean of a Pareto distribution is undefined if α ≤ 1 and grows arbitrarily large as α approaches 1 from above.
A more robust choice is the harmonic mean h B = l B (1 + 1/α), which is defined for all α > 0 (von . We use the harmonic mean in this article. We estimate α by fitting a Pareto distribution to the top two bins and calculating the maximum likelihood estimate from the following formula (Quandt, 1966): If the survey provides the grand mean µ, then we need not assume that top-bin incomes follow a Pareto distribution. Instead, we can calculatê which would be the mean of the top bin if the means of the lower bins were the midpoints m b . Thenμ B can serve as a pseudo-midpoint for the top bin.
When estimated in this way, it occasionally happens (e.g., in 4 percent of US counties) that the top bin's meanμ B is slightly less than its lower bound l B . This is infelicitous, but µ B can still be used, and the resulting Gini estimates are not necessarily bad. 1 The midpoint methods described in this section are implemented by the rpme command for Stata (Duan and von Hippel, 2015), where rpme stands for "robust Pareto midpoint estimator" (von Hippel and Powers, 2017). Except for mean-matching, the approach is also implemented in the binequality package for R (Scarpino et al., 2017).

Fitting parametric distributions
The weakness of the midpoint method is that it treats income as discrete. An alternative is to model income as a continuous variable X that fits some parametric CDF F (X|θ).
Here θ is a vector of parameters, which can be estimated by iteratively maximizing the log likelihood: While any parametric distribution can be considered, in practice it is hard to fit a distribution unless the number of parameters is small compared to the number of well-populated bins. Most investigators favor 2-, 3-, and 4-parameter distributions from the generalized beta family (McDonald and Ransom, 2008), which includes the following 10 distributions: the log normal, the log logistic, the Pareto (type 2), the gamma and generalized gamma, the beta 2 and generalized beta (type 2), the Dagum, the Singh-Maddala, and the Weibull.
A priori it is hard to know which distribution, if any, will fit well. The fit of a distribution can be tested by the following goodness-of-fit likelihood ratio statistic (von Hippel et al.,

2015)
: wherel is the maximized log likelihood. Under the null hypothesis that the fitted distribution is the true distribution of income, G 2 would follow a chi-square distribution with and B >0 is the number of bins with nonzero counts. We reject the fit of a distribution if G 2 has p < .05 in the null distribution.
In empirical data it is common to reject every distribution in the generalized beta family (Bandourian et al., 2002;von Hippel et al., 2015). In addition, some distributions may fail to converge, or may converge on parameter values that imply that the mean or variance of the distribution is undefined .
A solution is to fit all 10 distributions, screen out any with undefined moments and, among the distributions remaining, select the one that fits best according to the Akaike or 2 This formula includes the top bin B, where u B = ∞ and F (u B ) = 1. Some older articles (e.g., (McDonald and Ransom, 1979;McDonald and Xu, 1995;Bandourian et al., 2002)) add the following constant to the log likelihood: ln(T !) − B b=1 ln(n b !). This is not wrong, but it is unnecessary since adding a constant does not change the parameter values at which the log likelihood is maximized (Edwards, 1972).

Bayes information criterion (AIC or BIC):
Or, instead of selecting a single best-fitting distribution, one can average estimates across several candidate distributions weighted proportionately to a function of the AIC or BIC (specifically exp(−AIC/2) or exp(−BIC/2)) an approach known as model averaging (Burnham and Anderson, 2004). In general, model averaging yields better estimates than model selection, but when modeling binned incomes the advantage of model averaging is negligible .
Although it sounds broad-minded to fit 10 different distributions, there is limited diversity in the generalized beta family. All the distributions in the generalized beta family are unimodal and skewed to the right. Some distributions are quite similar (e.g., Dagum and generalized beta), and others rarely fit well (e.g., log normal, log logistic, Pareto).
So in practice the range of viable and contrasting distributions in the generalized beta family is small; you can fit just 3 well-chosen distributions (e.g., Dagum, gamma, and generalized gamma) and get estimates almost as good as those obtained from fitting all 10 .
We use the fitted distributios to estimate income statistics such as the mean, variance, Gini, or Theil. Sometimes the income statistic is a simple function of the distributional parameters θ, but other times the function is unknown, or hard to calculate. As a general solution, it is easier to calculate income statistics by applying numeric integration to appropriate functions of the fitted distribution (McDonald and Ransom, 2008).
When the grand mean is available, the distributional parameters could in theory be constrained to match the grand mean as well as approximate the bin counts. This would be difficult, though, since for distributions in the generalized beta family the mean is a complicated nonlinear function of the parameters. We have not attempted to constrain our parametric distributions to match a known mean.
The parametric approaches described in this section are implemented in the binequality package for R (Scarpino et al., 2017) and the mgbe command for Stata (Duan and von Hippel, 2015). Here mgbe stands for "multi-model generalized beta estimator" (von Hippel et al.,

Interpolated CDFs
Since parametric distributions may fit poorly, an alternative is to define a flexible nonparametric density that fits the bin counts exactly. Our nonparametric approach is called CDF interpolation.
To understand CDF interpolation, consider that binned data define B discrete points on the empirical cumulative distribution function (CDF). The empirical CDF for Nantucket is given in the last column of Table 1, which shows that 5% of Nantucket households make less than $10,000, 8% make less than $15,000, and so on. Formally, at an income of 0 the empirical CDF isF (0) = 0, and at each bin's upper bound u b the empirical CDF is the observed fraction of incomes that are less than u b -i.e.,F (u b ) = (n 1 + n 2 + · · · + n b )/T . 3 Now to estimate a continuous CDFF (x), x > 0, we just connect the dots. That is, we define a continuous nondecreasing function that interpolates between the B discrete points of the empirical CDF.
Then the estimated probability density function (PDF) is just the derivative of the interpolated CDFF (x). Note that the estimated PDF "preserves areas" (Morandi and Costantini, 1989)-i.e., preserves bin counts. That is, within in each bin (P (l b < x < u b ) =F (u b ) − F (l b )), the total of the estimated density is equal to the observed fraction of incomes that are in that bin (n b /T ).
The shape of the PDF depends on the function that interpolates the CDF.
• If the CDF is interpolated by line segments, then the CDF is polygonal, and the PDF is a step function that is discontinuous at the bin boundaries.
• If the CDF is interpolated more smoothly, by a continuously differentiable monotone cubic spline, then the PDF is piecewise quadratic -i.e., continuous at the bin boundaries and quadratic between them.
There remains a question of how to shape the CDF in the top bin, which typically has no upper bound. In our implementation, the CDF of the top bin can be rectangular, exponential, or Pareto. 4 Each of these distributions has one parameter, which is estimated as follows.
• If the grand mean of income is known, then the parameter shaping the top bin is constrained so that the mean of the fitted distribution matches the grand mean. It occasionally happens (e.g., in 4 percent of US counties) that we cannot reproduce the known mean this way, because the known mean is already less than the mean of the lower B − 1 bins without the tail. In that case, we make an ad hoc adjustment by shrinking the bin boundaries toward the origin -that is, by replacing (l b , u b ) with (sl b , su b ), where the shrinkage factor s < 1 is chosen so that a small tail can be added to reproduce the grand mean. The shrinkage factor is rarely less than .995.
• If the mean income is not known, we substitute an ad hoc estimate. We obtain that estimate by temporarily setting the upper bound of the top bin to u B = 2l B and calculating the mean of a step PDF fit to all B bins. Then we unbound the top bin and proceed as though the mean were known. Income statistics are estimated by applying numerical integration to functions of the fitted PDF or CDF.
The methods in this section are implemented by our binsmooth package for R (Hunter and Drown, 2016). Within the binsmooth package, the stepbins function implements a step-function PDF (and polygonal CDF), while the splinebins function implements a cubic spline CDF (and piecewise quadratic PDF).

Recursive subdivision
Another way to obtain a smooth PDF estimate that preserves bin areas is to subdivide the bins into smaller bins, and then adjust the heights of the subdivided bins to shorten the jumps at the bin boundaries. This method, recursive subdivision, is implemented by the recbin command in the binsmooth package for R (Hunter and Drown, 2016). Recursive subdivision is slower and more computationally intensive than CDF interpolation and the resulting estimates are practically identical. We present the details of recursive subdivision in the Appendix.

Data and Results
Between 2006  The Census also published means and Gini coefficients for each county (Bee, 2012).
These statistics were estimated from exact incomes before binning, and so are more accurate than any estimate that could be calculated from the binned data. They are sample estimates which may differ from population values, but they remain a useful standard of comparison for our binned-data estimates. Table 1 summarized the binned incomes for Nantucket County. Figure 1 fits several distributions to the Nantucket data.

Results for Nantucket
The midpoint method is illustrated by gray spikes at the bin midpoints; the spike heights are proportional to the bin counts. The black step function is the PDF implied by a linear interpolation of the CDF, and the blue curve is the piecewise quadratic PDF implied by a cubic spline of the CDF. Both CDF interpolations fit the bin counts perfectly; in fact, their jagged appearance suggests they may overfit the data-a concern that we will revisit in the Conclusion. The step PDF looks slightly less volatile than the piecewise quadratic PDF, suggesting that the step PDF may be less overfit. The Dagum distribution is drawn in purple, the quadratic spline in blue, and the step function in black.
Gray spikes illustrate the midpoint method.
The purple curve is the Dagum distribution. The Dagum fits Nantucket better than other distributions from the generalized beta family, but it does not fit well. It fails the G 2 goodness-of-fit test, and visually it fits the bin counts poorly. For example, between $70,000 and $150,000 the Dagum curve suggests there should be substantially fewer households than there are, and above $150,000 it suggests that there should be more.
Except for the Dagum distribution, all the methods in Figure 1 are calibrated to reproduce the grand mean. Table 2 summarizes the Nantucket estimates. The true mean is $137,000 and the true Gini is .547. All the methods underestimate these quantities. When fit without knowledge of the true mean, every method underestimates the mean by 12 to 20 percent and the Gini by 15 to 21 percent, with the simple midpoint method coming closer than its more sophisticated competitors.
When given the true mean, the midpoint and CDF interpolation methods do much better. They still underestimate the Gini, but only by 2 to 7 percent. The closest estimate is obtained by linear interpolation of the CDF. A smoother cubic spline interpolation does a little worse, but still better than the midpoint method. Although the estimates for Nantucket are less accurate overall than the estimates for most other counties, the relative performance of different methods in Nantucket is similar to what we will see elsewhere. We can summarize the accuracy of Gini estimates in several ways. For a single county j, the percent estimation error is e j = 100 × (θ j − θ j )/θ j . Then across all counties, the percent relative bias is the mean of e j , the percent root mean squared error (RMSE) is the square root of the mean of e 2 j , and the reliability is the squared correlation between θ j and θ j . Table 3 summarizes our findings. If the estimators ignore the published county means, then estimated Ginis have biases between 0% and -3%, RMSEs between 3% and 4%, and reliabilities between 82% and 88%. The interpolated CDF estimates have the best bias, the best RMSE, and the second best reliability, and they are just as good with linear Accuracy of Gini estimates interpolation as with cubic spline interpolation. The parametric estimates have the best reliability, but the worst bias, the worst RMSE, and by far the worst runtime, at 4.5 hours.

Results for all US counties
When the methods are constrained to match the published county means, the estimates improve dramatically. The bias shrinks to 0-1%, the RMSE shrinks to 1-2%, and the reliability grows to 98-99%. The midpoint estimates are excellent, and the interpolated CDF estimates are even better, and just as good with linear interpolation as with cubic spline interpolation.
The differences among the methods are much smaller than the improvement that comes from constraining any method to match the mean. Of course, this observation is only helpful when the mean is known. GHz.

Conclusion
CDF interpolation produces estimates that are at least a little better than midpoint or parametric estimates, whether the true mean is known or not. And CDF interpolation runs much faster than parametric estimation, thought not as fast as midpoint estimation.
We initially suspected that cubic spline interpolation would improve on simple linear interpolation, but empirically this turns out to be false. In estimating county Ginis, linear CDF interpolation was at least as accurate as cubic spline interpolation.
The accuracy of linear CDF interpolation is remarkable, since it implies a step function for the PDF.
Step PDFs look clearly unrealistic, especially in the top and bottom bins where the step function is flat while the true distribution likely has an upward or downward slope (Cloutier, 1995). Our step PDF permits a downward Pareto or exponential slope in the top bin, and our cubic spline CDF can fit an upward slope to the bottom bin. But neither of these refinements does much to improve the accuracy of Gini estimates.
The differences in accuracy among the methods are small, and they are dwarfed by the improvement in accuracy that comes from knowing the grand mean. By constraining binned-data methods to match a known mean, we can typically get county Gini estimates that are within 1-2% of the estimates we would get if the data were not binned. Our binsmooth package for R can constrain interpolated CDFs to match a known mean, and our rpme command for Stata can constrain the top-bin midpoint to match a known mean as well. We have not constrained our parametric distributions to match a known mean; we believe it would be difficult to do so.
While the mean-constrained estimates are very accurate, there may be room for improvement when the mean is unknown. Perhaps the most promising idea is additional smoothing. As we noticed in Figure fig:comparepdfs, interpolated CDFs can be a bit jagged and may "overfit" the sample in the sense that they find nooks and crannies that might not appear in another sample or in the population. Likewise interpolated CDFs may be overfit to a specific set of bin boundaries. If the fitted CDF were a little smoother and did not quite preserve the counts of the least populous bins, it might fit the population and other samples (perhaps with different bin boundaries) a little better.
A Appendix: PDF smoothing by recursive subdivision Recursive subdivision is another way to smooth the fitted PDF. Like CDF interpolation, recursive subdivision preserves bin areas, but recursive subdivision is more computationally intensive and produces very similar results. Recursive subdivision is implemented by the recbins function in our binsmooth package for R. A slight change of notation will be helpful. Since the upper bound u b of each bin is equal to the lower bound l b+1 of the next, we can think the bins as having a set of "edges" Start by fitting a step PDF. Let h b be the height of the step PDF in the bin [e b , e b+1 ).
Given parameters ε 1 ∈ (0, 0.5) and ε 2 ∈ (0, 1), the subdivision process begins by introducing new bin edges l and r between e b and e b+1 such that (l + r)/2 = (e b + e b+1 )/2 and r − l = In order for the above formulas to apply to the top and bottom bins, we create pseudobins above and below them with heights of zero. This ensures that the subdivided PDF will tend toward a height of zero at the lower edge of the bottom bin and the upper edge of the top bin.
The smoothed PDF is obtained from the step PDF by applying the subdivision process to each bin, then applying the process again to each subdivided bin, and so on, until the desired level of smoothness is reached. In practice, three rounds of subdivision are sufficient to produce a reasonably smooth PDF, and we found that choosing ε 1 = 0.25 and ε 2 = 0.75 Figure 3: Bin subdivision. Each original bin is replaced by three new bins (bold, dashed) such that the bin area is preserved.
produced nicely smoothed PDF's from most empirical data sets. Figure 4 shows the result of recursive subdivision in Nantucket.
Unfortunately, if the original step PDF was constrained to match a known mean, the subdivision process may cause the mean to deviate slightly. But the estimated Gini typically remains quite accurate.