Have the Nones Leveled Off?

June 29, 2024 AllenDowney

Last month Ryan Burge published “The Nones Have Hit a Ceiling“, using data from the 2023 Cooperative Election Study to show that the increase in the number of Americans with no religious affiliation has hit a plateau. Comparing the number of Atheists, Agnostics, and “Nothing in Particular” between 2020 and 2023, he found that “the share of non-religious Americans has stopped rising in any meaningful way.”

When I read that, I was frustrated that the HERI Freshman Survey had not published new data since 2019. I’ve been following the rise of the “Nones” in that dataset since one of my first blog articles.

As you might guess, the Freshman Survey reports data from incoming college students. Of course, college students are not a representative sample of the U.S. population, and as rates of college attendance have increased, they represent a different slice of the population over time. Nevertheless, surveying young adults over a long interval provides an early view of trends in the general population.

Well, I have good news! I got a notification today that HERI has published data tables for the 2020 through 2023 surveys. They are in PDF, so I had to do some manual data entry, but I have results!

Religious preference

Among other questions, the Freshman Survey asks students to select their “current religious preference” from a list of seventeen common religions, “Other religion,” “Atheist”, “Agnostic”, or “None.”

The options “Atheist” and “Agnostic” were added in 2015. For consistency over time, I compare the “Nones” from previous years with the sum of “None”, “Atheist” and “Agnostic” since 2015.

The following figure shows the fraction of Nones from 1969, when the question was added, to 2023, the most recent data available.

The blue line shows data until 2015; the orange line shows data from 2015 through 2019. The gray line shows a quadratic fit. The light gray region shows a 95% predictive interval.

The quadratic model continues to fit the data well and the recent trend is still increasing, but if you look at only the last few data points, there is some evidence that the rate of increase is slowing.

But not for women

Now here’s where things get interesting. Until recently, female students have been consistently more religious than male students. But that might be changing. The following figure shows the percentages of Nones for male and female students (with a missing point in 2018, when this breakdown was not available).

Since 2019, the percentage of Nones has increased for women and decreased for men, and it looks like women may now be less religious. So the apparent slowdown in the overall trend might be a mix of opposite trends in the two groups.

The following graph shows the gender gap over time, that is, the difference in percentages of male and female students with no religious affiliation.

The gap was essentially unchanged from 1990 to 2020. But in the last three years it has changed drastically. It now falls outside the predictive range based on past data, which suggests a change this large would be unlikely by chance.

Similarly with attendance at religious services, the gender gap has closed and possibly reversed.

UPDATE: Ryan Burge looked at the gender gap in CES and GSS data and found similar results: especially among young people, the gender gap has either disappeared or crossed over. And Ryan pointed me to this article by Dan Cox and Kelsey Eyre Hammond which reports similar trends in data from the Survey Center on American Life.

Attendance

The survey also asks students how often they “attended a religious service” in the last year. The choices are “Frequently,” “Occasionally,” and “Not at all.” Respondents are instructed to select “Occasionally” if they attended one or more times, so a wedding or a funeral would do it.

The following figure shows the fraction of students who reported any religious attendance in the last year, starting in 1968. I discarded a data point from 1966 that seems unlikely to be correct.

There is a clear dip in 2021, likely due to the pandemic, but the last two data points have returned to the long-term trend.

Data Source

The data reported here are available from the HERI publications page. Since I entered the data manually from PDF documents, it’s possible I have made errors.

Should divorce be more difficult?

June 14, 2024 AllenDowney

“The Christian right is coming for divorce next,” according to this recent Vox article, and “Some conservatives want to make it a lot harder to dissolve a marriage.”

As always when I read an article like this, I want to see data — and the General Social Survey has just the data I need. Since 1974, they have asked a representative sample of the U.S. population, “Should divorce in this country be easier or more difficult to obtain than it is now?” with the options to respond “Easier”, “More difficult”, or “Stay as is”.

Here’s how the responses have changed over time:

Since the 1990s, the percentage saying divorce should be more difficult has dropped from about 50% to about 30%. [The last data point, in 2022, may not be reliable. Due to disruptions during the COVID pandemic, the GSS changed some elements of their survey process — in the 2021 and 2022 data, responses to several questions have deviated from long-term trends in ways that might not reflect real changes in opinion.]

If we break down the results by political alignment, we can see whether these changes are driven by liberals, conservatives, or both.

Not surprisingly, conservatives are more likely than liberals to believe that divorce should be more difficult, by a margin of about 20 percentage points. But the percentages have declined in all groups — and fallen below 50% even among self-described conservatives.

As the Vox article documents, conservatives in several states have proposed legislation to make divorce more difficult. Based on the data, these proposals are likely to be unpopular.

To see my analysis, you can run this notebook on Colab. For similar analysis of other topics, see Chapter 11 of Probably Overthinking It.

Which Standard Deviation?

June 8, 2024 AllenDowney

It’s another installment in Data Q&A: Answering the real questions with Python. Previous installments are available from the Data Q&A landing page.

standard_dev

Which Standard Deviation?¶

Here’s a question from the Reddit statistics forum.

When do we use N and when N-1 for computation of sd and cov?

So I was doing a task, where I had a portfolio of shares and I was given their yields. Then I had to calculate [covariance], [standard deviation] etc. But I simply do not know when to use N and when N-1. I only know that it has to do something with degrees of freedom. Can someone explain it to me like to a 10 year old? Thanks!

If you look up the formula for standard deviation, you are likely to find two versions. One has the sample size, N, in the denominator; the other has N-1.

Sometimes the explanation of when you should use each of them is not clear. And to make it more confusing, some software packages compute the N version by default, and some the N-1 version.

Let’s see if we can straighten it all out.

Click here to run this notebook on Colab.

I’ll download a utilities module with some of my frequently-used functions, and then import the usual libraries.

In [1]:

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("Downloaded " + str(local))
    return filename

download('https://github.com/AllenDowney/DataQnA/raw/main/nb/utils.py')

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

from utils import decorate

In [2]:

# install the empiricaldist library, if necessary

try:
    import empiricaldist
except ImportError:
    !pip install empiricaldist

Samples and Populations¶

Suppose we are given a sequence of numbers and asked to compute the standard deviation. As an example, we’ll use the sequence from 0 to 9.

In [3]:

N = 10
data = np.arange(N)
data

Out[3]:

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

The first step is to compute the mean.

In [4]:

m = np.mean(data)
m

Out[4]:

4.5

The deviations are the distances of each point from the mean.

In [5]:

deviations = data - m
deviations

Out[5]:

array([-4.5, -3.5, -2.5, -1.5, -0.5,  0.5,  1.5,  2.5,  3.5,  4.5])

Next we compute the sum of the squared deviations.

In [6]:

ssd = np.sum(deviations **2)
ssd

Out[6]:

82.5

Then we compute the variance — we’ll start with the version that has N in the denominator.

In [7]:

var = ssd / N

Finally, the standard deviation is the square root of variance.

In [8]:

std = np.sqrt(var)
std

Out[8]:

2.8722813232690143

And here’s the version with N-1 in the denominator.

In [9]:

var = ssd / (N-1)
std = np.sqrt(var)
std

Out[9]:

3.0276503540974917

With N=10, the difference between the version is non-negligible, but this is pretty much the worst case. With larger sample sizes, the difference in smaller — and with smaller sample sizes, you probably shouldn’t compute a standard deviation at all.

By default, NumPy computes the N version.

In [10]:

np.std(data)

Out[10]:

2.8722813232690143

But with the optional argument ddof=1, it computes the N-1 version.

In [11]:

np.std(data, ddof=1)

Out[11]:

3.0276503540974917

By default, Pandas computes the N-1 version.

In [12]:

pd.Series(data).std()

Out[12]:

3.0276503540974917

But with the optional argument ddof=0, it computes the N version.

In [13]:

pd.Series(data).std(ddof=0)

Out[13]:

2.8722813232690143

It is not ideal that the two libraries have different default behavior. And it might not be obvious why the parameter that controls this behavior is called ddof.

The answer is related to OP’s question about “degrees of freedom”. To understand that term, suppose I ask you to think of three numbers. You are free to choose any first number, any second number, and any third number. In that case, you have three degrees of freedom.

Now suppose I ask you to think of three numbers, but they are required to add up to 10. You are free to choose the first and second numbers, but then the third number is determined by the requirement. So you have only two degrees of freedom.

If we are given a dataset with N elements, we generally assume that it has N degrees of freedom unless we are told otherwise. That’s what the ddof parameter does. It stands for “delta degrees of freedom”, where “delta” indicates a change or a difference. In this case, it is the difference between the presumed degrees of freedom, N, and the degrees of freedom that should be used for the computation, N - ddof.

So, with ddof=0, the denominator is N. With ddof=1, the denominator is N-1.

Which one is right?¶

So when should you use one or the other? It depends on whether you are describing data or making a statistical inference.

The N version is a descriptive statistic — it is quantifies the variability of the data.
The N-1 version is an estimator — if data is a sample from a population, we can use it to infer the standard deviation of the population.

To be more specific, the N-1 version is an almost unbiased estimator, which means that it gets the answer almost right, on average, in the long run. Let’s see how that works.

Consider a uniform distribution from 0 to 10.

In [14]:

from scipy.stats import uniform

dist = uniform(0, 10)

Here are the actual mean and standard deviation of this distribution, computed analytically.

In [15]:

dist.mean(), dist.std()

Out[15]:

(5.0, 2.8867513459481287)

If we generate a large random sample from this distribution, we expect its standard deviation to be close to 2.9.

In [16]:

sample = dist.rvs(size=10000)
sample.std(ddof=0), sample.std(ddof=1)

Out[16]:

(2.8816683774305085, 2.881812471656537)

It is, and with a large sample size, the difference between the N version and the N-1 version is negligible.

But let’s see what happens if we generate small samples many times. First we’ll compute the N version of the standard deviation for 10001 samples.

In [17]:

standard_deviations_N = [dist.rvs(N).std(ddof=0) for i in range(10001)]

Ideally, the mean of these estimates should be close to the actual mean of the distribution, which is about 2.9.

In [18]:

np.mean(standard_deviations_N)

Out[18]:

2.6966272770954323

But it’s not — it is too low on average, which means it is a biased estimator.

Let’s see if the N-1 version is better.

In [19]:

standard_deviations_Nminus1 = [dist.rvs(N).std(ddof=1) for i in range(10001)]

np.mean(standard_deviations_Nminus1)

Out[19]:

2.8482429099753923

Yes! The N-1 version is a much less biased estimator of the standard deviation.

That’s why the N-1 version is called the “sample standard deviation”, because it is appropriate when we are using a sample to estimate the standard deviation of a population.

If, instead, we are able to measure an entire population, we can use the N version — which is why it is called the “population standard deviation”.

Should we care?¶

On one hand, it is impressive that such a simple correction yields a much better estimator. On the other hand, it almost never matters in practice.

If you have a large sample, the difference between the two versions is negligible.
If you have a small sample, you can’t make a precise estimate of the standard deviation anyway.

To demonstrate the second point, let’s look at the distributions of the estimates using the two versions.

In [20]:

sns.kdeplot(standard_deviations_N, label='N version')
sns.kdeplot(standard_deviations_Nminus1, label='N-1 version')

decorate(xlabel='Estimated standard deviation')

No description has been provided for this image

Compared to the variation in the estimates, the difference between the versions is small.

In summary, if your sample size is small and it is really important to avoid underestimating the standard deviation, use the N-1 correction. Otherwise it doesn’t matter.

In [23]:

data = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
np.std(data)

Out[23]:

2.8722813232690143

In [24]:

pd.Series(data).std()

Out[24]:

3.0276503540974917

Data Q&A: Answering the real questions with Python

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

What is a percentile rank?

June 1, 2024 AllenDowney

It’s another installment in Data Q&A: Answering the real questions with Python. Previous installments are available from the Data Q&A landing page.

percentile_rank

What is a Percentile Rank?¶

Here’s a question from the Reddit statistics forum.

What’s the difference between cumulative frequency as a percentage and percentiles?

Example that got me super confused:

[Here’s the data provided by OP]:

In [1]:

import pandas as pd

columns = ['score', 'frequency', 'cumulative frequency', 'cumulative percentage']

data = [
    (8, 1, 1, 2),
    (9, 1, 2, 4),
    (10, 1, 3, 6),
    (11, 3, 6, 13),
    (12, 5, 11, 23),
    (13, 7, 18, 38),
    (14, 8, 26, 54),
    (15, 9, 35, 73),
    (16, 6, 41, 85),
    (17, 4, 45, 94),
    (18, 2, 47, 98),
    (19, 1, 48, 100),
]
df = pd.DataFrame(data, columns=columns)
df

Out[1]:

	score	frequency	cumulative frequency	cumulative percentage
0	8	1	1	2
1	9	1	2	4
2	10	1	3	6
3	11	3	6	13
4	12	5	11	23
5	13	7	18	38
6	14	8	26	54
7	15	9	35	73
8	16	6	41	85
9	17	4	45	94
10	18	2	47	98
11	19	1	48	100

OP continues:

To me, the last [column], the “cumulative frequency” expressed in percentages is a percentile: it’s the percentage of scores lower then or the same as the score it’s calculated for.

However, in my class it’s explained that the percentile for score 14 would be 46, and that the percentile for score 17 would be 90, not 94.

The way they arrive at this is that upper and lower limits for the percentile calculation have to be set, in the case of score 14 that would be 14 and 13, and then the percentile would be calculated as (upper limit+lower limit)/2.

In the case of 14: (38+54)/2 = 46

This makes absolutely no sense to me: why are we introducing arbitrary limits to average out when the cumulative frequency in percentages, to me, seems to meet the definition of a percentile perfectly?

What’s at issue here is the definition of percentile and percentile rank. As an example, let’s consider the distribution of height, and suppose the median is 168 cm. In that case, 168 cm is the 50th percentile, and if someone is 168 cm tall, their percentile rank is 50%.

By this definition, the percentile rank for a particular quantity is its cumulative frequency expressed as a percentage, exactly as OP suggests. If you are taller than, or the same height as, 50% of the population, your percentile rank is 50%.

However, some classes teach the alternative definition OP presents. So, let me explain the two definitions, and we can discuss the pros and cons. At the risk of giving away the ending, here is my conclusion:

Percentiles and percentile ranks are perfectly well defined, and I don’t think we should encourage variations on the definitions.
In a dataset with many repeated values, percentile ranks might not measure what we want — in that case we might want to compute a different statistic, but then we should give it a different name.

In other words, I agree with OP.

Click here to run this notebook on Colab.

I’ll download a utilities module with some of my frequently-used functions, and then import the usual libraries.

In [2]:

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("Downloaded " + str(local))
    return filename

download('https://github.com/AllenDowney/DataQnA/raw/main/nb/utils.py')

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

from utils import decorate

In [3]:

# install the empiricaldist library, if necessary

try:
    import empiricaldist
except ImportError:
    !pip install empiricaldist

Percentiles and percentile ranks¶

From the dataset OP provided, we can extract the scores and their frequencies.

In [4]:

pairs = df[['frequency', 'score']].values

And we can reconstitute the sample by making the given number of copies of each score.

In [5]:

sample = np.concatenate([[score] * freq for freq, score in pairs])
sample

Out[5]:

array([ 8,  9, 10, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13,
       13, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15,
       15, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 19])

Now suppose we want to compute the median score and the interquartile range (IQR). We can use the NumPy function percentile to compute the 25th, 50th, and 75th percentiles.

In [6]:

np.percentile(sample, [25, 50, 75])

Out[6]:

array([13., 14., 16.])

So the median score is 14 and the IQR is 16 - 13 = 3.

Going in the other direction, if we are given a score, q, we can compute its percentile rank by computing the fraction of scores less than or equal to q, expressed as a percentage.

In [7]:

q = 14
np.mean(sample <= q) * 100

Out[7]:

54.166666666666664

If someone gets the median score of 14, their percentile rank is about 54%.

Here we see the first thing that bothers people about percentiles and percentile ranks: with a finite dataset, they are not invertible. If you compute the percentile rank of the 50th percentile, the answer is not always 50%. To see why, let’s consider the CDF.

The CDF¶

Percentiles and percentile ranks are closely related to the the cumulative distribution function (CDF). To demonstrate, we can use empiricaldist to make a Cdf object from the reconstituted sample.

In [8]:

from empiricaldist import Cdf

cdf = Cdf.from_seq(sample)
cdf

Out[8]:

	probs
8	0.020833
9	0.041667
10	0.062500
11	0.125000
12	0.229167
13	0.375000
14	0.541667
15	0.729167
16	0.854167
17	0.937500
18	0.979167
19	1.000000

A Cdf object is a Pandas Series that contains the observed quantities as an index and their cumulative probabilities as values. We can use square brackets to look up a quantity and get a cumulative probability.

In [9]:

cdf[14]

Out[9]:

0.5416666666666666

If we multiply by 100, we get the percentile rank.

But square brackets only work with quantities in the dataset. If we look up any other quantity, that’s an error. However, we can use parentheses to call the Cdf object like a function.

In [10]:

cdf(14)

Out[10]:

array(0.54166667)

And that works with any numerical quantity.

In [11]:

cdf(13.5)

Out[11]:

array(0.375)

Cdf provides a step method that plots the CDF as a step function, which is what it technically is.

In [12]:

cdf.step(label='')
decorate(ylabel='CDF')

To be more explicit, we can put markers at the top of each vertical segment to indicate how the function is evaluated at one of the observed quantities.

In [13]:

cdf.step(label='')
cdf.plot(style='o', label='')
decorate(ylabel='CDF')

The Cdf object provides a method called inverse that computes the inverse CDF — that is, if you give it a cumulative probability, it computes the corresponding quantity.

In [14]:

p1 = 0.375
cdf.inverse(p1)

Out[14]:

array(13.)

The result from an inverse lookup is sometimes called a “quantile”, which is the name of a Pandas method that computes quantiles.

If we put the sample into a Series, we can use quantile to compute the quantity to look up the cumulative probability p1.

In [15]:

sample_series = pd.Series(sample)
sample_series.quantile(p1)

Out[15]:

13.625

By default, quantile uses interpolates linearly between the observed values. If we want this function to behave according to the definition of the CDF, we have to specify a different kind of interpolation.

In [16]:

sample_series.quantile(p1, interpolation='lower')

Out[16]:

If the result from the inverse CDF is called a quantile, you might wonder what we call the result from the CDF. By analogy with percentile and percentile rank, I think it should be called a quantile rank, but no one calls it that. As far as I can tell, it’s just called a cumulative probability.

So what’s wrong with percentile rank?¶

Percentiles and percentile ranks have a perfectly good definition, which follows from the definition of the CDF. So what’s the problem? Well, I have to admit — the definition is a little arbitrary.

In the example, suppose your score is 14. Your score is strictly better than 37.5% of the other scores.

In [17]:

less = np.mean(sample < 14)
less

Out[17]:

0.375

It’s strictly worse than 45.8%

In [18]:

more = np.mean(sample > 14)
more

Out[18]:

0.4583333333333333

And equal to 16.7%.

In [19]:

same = np.mean(sample == 14)
same

Out[19]:

0.16666666666666666

So if we want a single number that quantifies your performance relative to the rest of the class, which number should we use?

The definition of the CDF suggests we should report the fraction of the class whose score is less than or equal to yours. But that is an arbitrary choice. We could just as easily report the fraction whose score is strictly less — or the midpoint of these extremes. For a score of 14, here’s the midpoint:

In [20]:

less + same / 2

Out[20]:

0.4583333333333333

That’s the result OP’s teacher was expecting, and to be fair, that’s Wikipedia’s definition of percentile rank. But I don’t like it.

Discussion¶

I prefer the CDF-based definition of percentile rank because it’s consistent with the way most computational tools work. The midpoint-based definition feels like a holdover from the days of working with small datasets by hand.

That’s just my preference — if people want to compute midpoints, I won’t stop them. But for the sake of clarity, we should give different names to different statistics.

Historically, I think the CDF-based definition has the stronger claim on the term “percentile rank”. For the midpoint-based definition, ChatGPT suggests “midpoint percentile rank” or “average percentile rank”. Those seem fine to me, but it doesn’t look like they are widely used.

Data Q&A: Answering the real questions with Python

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Probably Overthinking It

Data science, Bayesian Statistics, and other ideas

Browsed by
Month: June 2024

Have the Nones Leveled Off?

June 29, 2024 AllenDowney

Religious preference

But not for women

Attendance

Data Source

Should divorce be more difficult?

June 14, 2024 AllenDowney

Which Standard Deviation?

June 8, 2024 AllenDowney

Which Standard Deviation?¶

Samples and Populations¶

Which one is right?¶

Should we care?¶

What is a percentile rank?

June 1, 2024 AllenDowney

What is a Percentile Rank?¶

Percentiles and percentile ranks¶

The CDF¶

So what’s wrong with percentile rank?¶

Discussion¶

	score	frequency	cumulative frequency	cumulative percentage
0	8	1	1	2
1	9	1	2	4
2	10	1	3	6
3	11	3	6	13
4	12	5	11	23
5	13	7	18	38
6	14	8	26	54
7	15	9	35	73
8	16	6	41	85
9	17	4	45	94
10	18	2	47	98
11	19	1	48	100

	score	frequency	cumulative frequency	cumulative percentage
0	8	1	1	2
1	9	1	2	4
2	10	1	3	6
3	11	3	6	13
4	12	5	11	23
5	13	7	18	38
6	14	8	26	54
7	15	9	35	73
8	16	6	41	85
9	17	4	45	94
10	18	2	47	98
11	19	1	48	100

	score	frequency	cumulative frequency	cumulative percentage
0	8	1	1	2
1	9	1	2	4
2	10	1	3	6
3	11	3	6	13
4	12	5	11	23
5	13	7	18	38
6	14	8	26	54
7	15	9	35	73
8	16	6	41	85
9	17	4	45	94
10	18	2	47	98
11	19	1	48	100