Histograms: More Than Meets The Plot

The way I’ve taught histograms in the past is simple. You take a numeric variable, then you use ggplot2, and slap it onto your x aesthetic. Then, run your code chunk. Once you’ve done so, choose appropriate binwidths that show the relationship of the data frequencies. Then, you’d have your acceptable histogram. Instead of calling it a day right there, let’s start thinking about how we can define histograms in math notation.

1 The Anatomy of a Histogram

Each of the rectangles on a histogram is a bin. Bins are intervals defined as

\[ B_1 = [0, \frac{1}{m}), B_2 = [\frac{1}{m}, \frac{2}{m}), ..., B_m = [\frac{m-1}{m}, \frac{m}{m}). \]

Notice that the definition of bins in this way puts us onto the interval \([0,1]\). We will be assuming our draws are \(iid \sim f \in [0,1]\). Our binwidth parameter based on the above is \(h=\frac{1}{m}\). As you’ve probably noticed before, ggplot2 has these options in geom_histogram() for you to specify. Let’s start our study of histograms by taking a look at a histogram of max heart rates achieved by cardiology patients. In this dataset, our \(n=180\).

# * DECIDING MY BINWIDTH BASED ON VISUALS
our_range <- range(heart_data$resting_blood_pressure)[2] - range(heart_data$resting_blood_pressure)[1]
# * PLOTTING
ggplot(heart_data, aes(x=resting_blood_pressure)) +
  geom_histogram(binwidth=our_range/22, col="pink3", fill="pink") +
  xlab("resting blood pressure") +
  ggtitle("Max Heart Rates of Cardiology Patients") +
  theme_minimal()

When I teach, this is where I’d call it the “end of the day” for the histogram. We got to see an overall density pattern, and that’s all it’d be used for: visual purposes. In this blog, we’ll see how histograms are estimators for the true probability density of a random variable.

In the above histogram, \(h \approx 3.909\). This doesn’t seem like it quite fits the bill of \(h=\frac{1}{m}\), but this is because our \(X_1, ..., X_n\) are not \(iid \sim \, f \in [0,1]\). Let’s see what happens if we fix our data to be on the interval \([0,1]\).

# * RESCALING WITH MUTATE
heart_data <- heart_data %>% mutate(heart_rescale=max_heart_rate_achieved/max(max_heart_rate_achieved))
# * DECIDING MY NEW BINWIDTH BASED ON VISUALS
our_range_2 <- our_range / max(heart_data$resting_blood_pressure)
# * PLOTTING
ggplot(heart_data, aes(x=heart_rescale)) +
  geom_histogram(binwidth=our_range_2/22, col="pink3", fill="pink") +
  xlab("resting blood pressure") +
  ggtitle("Max Heart Rates of Cardiology Patients (Scaled)") +
  theme_minimal()

Having done the transformation, we are intrigued to see our histogram shape shift. Now, we can start diving into histogram estimators.

2 The Histogram Estimator

Let \(c_j\) be the number of observations in \(B_j\). We’ll let the probability of being in each bin be \(\hat{p}_j=c_j/n\). Finally, also let \(p_j=\int_{B_j} f(u) du\). Our histogram estimator is defined as follows

\[ \hat{f}_n(x)=\sum_{i=1}^{n} \frac{\hat{p}_j}{h} \textbf{1}_{\{x \in B_j\}}. \]

The histogram estimator \(\hat{f}_n(x)\) has expectation \(E(\hat{f}_n(x)) \approx f(x)\) where \(f(x)\) is the true probability density.

3 Expectation and Variance of the Histogram Estimator

The expectation of the histogram estimator is

\[ E(\hat{f}_n(x)) = \frac{p_j}{h} \] and

\[ \text{Var}(\hat{f}_n(x)) = \frac{p_j(1-p_j)}{nh^2} \]

where \(h\) is binwidth and \(n\) is the total number of bins.

4 Choosing Binwidth Based on the Bias-Variance Tradeoff

Recall that

\[ \text{Risk}(\theta)=\text{Bias}^2(\theta) + \text{Var}(\theta) \]

for some arbitrary estimator \(\theta\). Making use the expectation and variance of the histogram estimator, the bias-variance tradeoff, and the fact that for any \(u \in B_j\), \(f(u) \approx f(x)+(u-x)f'(x)\), we can calculate an optimal binwidth.

4.1 By Theorem

The following theorem shows that you can choose an optimal binwidth for your histogram estimator to minimize risk (mean integrated square error - MISE). The theorem shows that there does indeed exist an optimal binwidth, but it is based on the true distribution of \(f(u)\) which we don’t have in practice. This is Theorem 20.4 in Wasserman.

Given that \(\int [f'(u)]^2 du < \infty\), we have that

\[ \text{Risk}(\hat{f}_n, f) \approx \frac{h^2}{12} \int [f'(u)]^2 du + \frac{1}{nh}. \]

The value that of \(h\) (the binwidth size) that minimizes the above risk is

\[ h^*= \frac{1}{n^{1/3}} \bigg( \frac{6}{\int [f'(u)]^2 du} \bigg)^{1/3}. \]

When that binwidth is chosen, \(\text{Risk}(\hat{f}_n, f) \approx \frac{C}{n^{2/3}}\) where \(C=(3/4)^{2/3}\big( \int [f'(u)]^2 du \big)^{1/3}\). The MISE converges to 0 at a rate \(n^{-2/3}\).

4.2 By Estimation

In practice, we would be interested in estimating the risk function so that it doesn’t depend on \(f(u)\). The cross-validation estimator of risk (estimated risk) is

\[ \hat{J}(h) = \int [f'_n(x)]^2 dx - \frac{2}{n}\sum_{i=1}^{n} \hat{f}_{(-i)}(X_i). \] The \((-i)\) stands for the removal of the \(i^{th}\) observation. Also,

\[ \hat{J}(h) = \frac{2}{(n-1)h} - \frac{(n+1)}{(n-1)} \sum_{j=1}^{m} \hat{p}_j^2. \]

4.3 Estimating Risk for Heart Disease Data

Let’s look back at the heart disease data. To estimate risk, I am going to use Wasserman’s cv.hist.fun from the All of Statistics website. Source code is available here. The code implements the second formula for cross-validated risk estimate: \(\hat{J}(h) = \frac{2}{(n-1)h} - \frac{(n+1)}{(n-1)} \sum_{j=1}^{m} \hat{p}_j^2\).

# * USING WASSERMAN'S FUNCTION
source("histograms-1-source.R")
function_output <- cv.hist.fun(heart_data %>% pull(heart_rescale))
risk_estimates <- data.frame(cbind(h=function_output$h, n_bins=function_output$nbins, est_risk=function_output$risk))
risk_estimates %>% pull(est_risk) %>% summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.040  -2.470  -2.327  -2.333  -2.159  -1.794 

According to this output, having 180 bins minimizes estimated MISE. Since this matches the value of \(n\), it may lead you to think each bin is an observation. however, this definitely does not correspond to each observation having a bin. Recall that having 180 bins means

\[ B_1 = [0.4752, 0.4752 + \frac{1}{180}), B_2 = [0.4808, 0.4919), ..., B_m = [\frac{m-1}{m}, \frac{m}{m}). \]

where \(m=180\). This does not guarantee that there will only be one value in each of these bins.

h <- function_output$hbest
# * PLOTTING
ggplot(heart_data, aes(x=heart_rescale)) +
  geom_histogram(binwidth=h, col="pink3", fill="pink") +
  xlab("resting blood pressure") +
  ggtitle("Max Heart Rates of Cardiology Patients (Scaled)") +
  theme_minimal()

5 Confidence Bands

According to Wasserman’s 20.9 Definition, a pair of functions \(l(x)\) (lower) and \(u(x)\) (upper) is a \(1-\alpha\) confidence band if P(\(\bar{f}_n(x)\) of being within the two functions) \(\geq 1-\alpha\). That is,

\[ P\big(l(x) < \bar{f}_n(x) < u(x)\big) \geq 1-\alpha \]

where \(\bar{f}_n(x)=E(\hat{f}_n(x)) = \frac{p_j}{h} \, \forall x\in B_j\) and where \(p_j=\int_{B_j}f(u)du\).

Using a function from Wasserman, we calculate a 95% confidence band. The lines below stand for the “top” of the histogram. Imagine the diagram to be shaded below those lines. The middle line stands for \(\bar{f}_n(x)\).

# * USING WASSERMAN'S FUNCTION
source("histograms-1-source.R")
function_output_2 <- ci.fun(heart_data %>% pull(heart_rescale))
envelope_df <- data.frame(cbind(grid=function_output_2$Grid, u=function_output_2$U, l=function_output_2$L, f=function_output_2$F))
# * PLOTTING ENVELOPE
ggplot(envelope_df, aes(x=grid)) + 
  geom_line(stat="identity", aes(x=grid, y=u), col="pink3", lwd=1, alpha=0.6) +
  geom_line(stat="identity", aes(x=grid, y=f), col="pink2", lwd=2) +
  geom_line(stat="identity", aes(x=grid, y=l), col="pink", lwd=1, alpha=0.6) +
  theme_minimal() +
  xlab("heart rate maximums") +
  ylab("frequencies") +
  ggtitle("95% Confidence Band for Cardiology Patients' Heart Rate Maximums")

At least to me, histograms have a lot more life to them than they used to. If you had to remember one thing from this blog, it’s that a histogram can have an optimal number of bins that will minimize MISE. This is not to say that when you’re visualizing you shouldn’t choose bins like I did in my first two figures. Choose bins in this way when you are estimating a probability density.

6 References

Wasserman, L. (2013). All of Statistics. Springer Science & Business Media.

• Risk estimation code.

• Confidence band code.

7 Data Source

Download link here.

1. Hungarian Institute of Cardiology. Budapest: Andras Janosi, M.D.
   
2. University Hospital, Zurich, Switzerland: William Steinbrunn, M.D.
   
3. University Hospital, Basel, Switzerland: Matthias Pfisterer, M.D.
   
4. V.A. Medical Center, Long Beach and Cleveland Clinic Foundation: Robert Detrano, M.D., Ph.D.