Section 4.4 Percentile-Based Statistics
In Chapter 3 we computed the arithmetic mean, which identifies a central point in a distribution, and the standard deviation, which quantifies how spread out the distribution is. And in a previous exercise we computed skewness, which indicates whether a distribution is skewed left or right. One drawback of all of these statistics is that they are sensitive to outliers. A single extreme value in a dataset can have a large effect on mean, standard deviation, and skewness.
An alternative is to use statistics that are based on percentiles of the distribution, which tend to be more robust, which means that they are less sensitive to outliers. To demonstrate, let’s load the NSFG data again without doing any data cleaning.
from nsfg import read_stata
dct_file = "2002FemPreg.dct"
dat_file = "2002FemPreg.dat.gz"
preg = read_stata(dct_file, dat_file)
Recall that birth weight is recorded in two columns, one for the pounds and one for the ounces.
birthwgt_lb = preg["birthwgt_lb"]
birthwgt_oz = preg["birthwgt_oz"]
If we make a
Hist object with the values from birthwgt_oz, we can see that they include the special values 97, 98, and 99, which indicate missing data.
$ from empiricaldist import Hist
Hist.from_seq(birthwgt_oz).tail(5)
birthwgt_oz
14.0 475
15.0 378
97.0 1
98.0 1
99.0 46
Name: , dtype: int64
The
birthwgt_lb column includes the same special values; it also includes the value 51, which has to be a mistake.
$ Hist.from_seq(birthwgt_lb).tail(5)
birthwgt_lb
15.0 1
51.0 1
97.0 1
98.0 1
99.0 57
Name: , dtype: int64
Now let’s imagine two scenarios. In one scenario, we clean these variables by replacing missing and invalid values with
nan, and then compute total weight in pounds. Dividing birthwgt_oz_clean by 16 converts it to pounds in decimal.
birthwgt_lb_clean = birthwgt_lb.replace([51, 97, 98, 99], np.nan)
birthwgt_oz_clean = birthwgt_oz.replace([97, 98, 99], np.nan)
total_weight_clean = birthwgt_lb_clean + birthwgt_oz_clean / 16
In the other scenario, we neglect to clean the data and accidentally compute the total weight with these bogus values.
total_weight_bogus = birthwgt_lb + birthwgt_oz / 16
The bogus dataset contains only 49 bogus values, which is about 0.5% of the data.
$ count1, count2 = total_weight_bogus.count(), total_weight_clean.count()
diff = count1 - count2
diff, diff / count2 * 100
(np.int64(49), np.float64(0.5421553441026776))
Now let’s compute the mean of the data in both scenarios.
$ mean1, mean2 = total_weight_bogus.mean(), total_weight_clean.mean()
mean1, mean2
(np.float64(7.319680587652691), np.float64(7.265628457623368))
The bogus values have a moderate effect on the mean. If we take the mean of the cleaned data to be correct, the mean of the bogus data is off by less than 1%.
$ (mean1 - mean2) / mean2 * 100
np.float64(0.74394294099376)
An error like that might go undetected — but now let’s see what happens to the standard deviations.
$ std1, std2 = total_weight_bogus.std(), total_weight_clean.std()
std1, std2
(np.float64(2.096001779161835), np.float64(1.4082934455690173))
$ (std1 - std2) / std2 * 100
np.float64(48.83274403900607)
The standard deviation of the bogus data is off by almost 50%, so that’s more noticeable. Finally, here’s the skewness of the two datasets.
def skewness(seq):
"""Compute the skewness of a sequence
seq: sequence of numbers
returns: float skewness
"""
deviations = seq - seq.mean()
return np.mean(deviations**3) / seq.std(ddof=0) ** 3
$ skew1, skew2 = skewness(total_weight_bogus), skewness(total_weight_clean)
skew1, skew2
(np.float64(22.251846195422484), np.float64(-0.5895062687577697))
$ # how much is skew1 off by?
(skew1 - skew2) / skew2
np.float64(-38.74658112171128)
The skewness of the bogus dataset is off by a factor of almost 40, and it has the wrong sign! With the outliers added to the data, the distribution is strongly skewed to the right, as indicated by large positive skewness. But the distribution of the valid data is slightly skewed to the left, as indicated by small negative skewness.
These results show that a small number of outliers have a moderate effect on the mean, a strong effect on the standard deviation, and a disastrous effect on skewness.
An alternative is to use statistics based on percentiles. Specifically:
-
The median, which is the 50th percentile, identifies a central point in a distribution, like the mean.
-
The interquartile range, which is the difference between the 25th and 75th percentiles, quantifies the spread of the distribution, like the standard deviation.
-
The quartile skewness uses the quartiles of the distribution (25th, 50th, and 75th percentiles) to quantify the skewness.
The
Cdf object provides an efficient way to compute these percentile-based statistics. To demonstrate, let’s make a Cdf object from the bogus and clean datasets.
cdf_total_weight_bogus = Cdf.from_seq(total_weight_bogus)
cdf_total_weight_clean = Cdf.from_seq(total_weight_clean)
The following function takes a
Cdf and uses its inverse method to compute the 50th percentile, which is the median (at least, it is one way to define the median of a dataset).
def median(cdf):
m = cdf.inverse(0.5)
return m
Now we can compute the median of both datasets.
$ median(cdf_total_weight_bogus), median(cdf_total_weight_clean)
(array(7.375), array(7.375))
The results are identical, so in this case, the outliers have no effect on the median at all. In general, outliers have a smaller effect on the median than on the mean.
The interquartile range (IQR) is the difference between the 75th and 25th percentiles. The following function takes a
Cdf and returns the IQR.
def iqr(cdf):
low, high = cdf.inverse([0.25, 0.75])
return high - low
And here are the interquartile ranges of the two datasets.
$ iqr(cdf_total_weight_bogus), iqr(cdf_total_weight_clean)
(np.float64(1.625), np.float64(1.625))
In general, outliers have less effect on the IQR than on the standard deviation — in this case they have no effect at all.
Finally, here’s a function that computes quartile skewness, which depends on three statistics:
-
The median,
-
The midpoint of 25th and 75th percentiles, and
-
The semi-IQR, which is half of the IQR.
def quartile_skewness(cdf):
low, median, high = cdf.inverse([0.25, 0.5, 0.75])
midpoint = (high + low) / 2
semi_iqr = (high - low) / 2
return (midpoint - median) / semi_iqr
And here’s the quartile skewness for the two datasets.
$ qskew1 = quartile_skewness(cdf_total_weight_bogus)
qskew2 = quartile_skewness(cdf_total_weight_clean)
qskew1, qskew2
(np.float64(-0.07692307692307693), np.float64(-0.07692307692307693))
The small number of outliers in these examples has no effect on the quartile skewness. These examples show that percentile-based statistics are less sensitive to outliers and errors in the data.
