Skip to content
10 min read csv pyarrow grouping stack-unstack pivot-table

Bamboo Weekly #55: IVF (solutions)

Get better at: CSV, PyArrow, grouping, stack-unstack, and pivot tables.

Bamboo Weekly #55: IVF (solutions)

Two administrative notes:

  1. I’ll be holding office hours for paid subscribers this coming Sunday, March 3rd. Look for a notice tomorrow with the full Zoom info.
  2. I’m giving two talks at PyCon US in May, including one about dates and times in Pandas. I’m also a sponsor of the conference, and will have a booth there telling people about my course subscription program, corporate training, and (of course) Bamboo Weekly. I hope to see you there!

This week, we looked at in vitro fertilization (IVF). The widespread fertility treatment made headlines when Alabama’s Supreme Court ruled that embryos are children, as per their state constitution. IVF typically involves the fertilization of several embryos outside of the womb; those that aren’t used are then frozen or destroyed. As a result, the court’s ruling has effectively shut down IVF treatment in Alabama. It has also sent political shock waves through much of the US, and might play a role in US elections this coming November.

This week’s data set came from the Centers for Disease Control and Prevention (CDC). The CDC surveyed a large number of fertility clinics across the US about their use of IVF. The questions ranged from the number of retrievals (i.e., of egg cells) that were performed to the number of transfer (i.e., insertion of fertilized embryos) that were done, to the transfer technique used, to the number of births (single and multiple) that resulted from IVF. The data set was from October 2023, with data reported in 2021; from what I can tell, IVF research commonly lags by about two years, because of the time involved in achieving pregnancy and birth.

Data and six questions

The home page for this data set is here:

https://data.cdc.gov/Assisted-Reproductive-Technology-ART-/2021-Final-Assisted-Reproductive-Technology-ART-Su/9tjt-seye/about_data

Note that while the above URL is where you have to get the data, it took me a while to figure out how to get it. Here's what I did:

  1. Click on the "export" button at the top of the page
  2. Select "CSV file" for export format (which should be selected anyway, as the default)
  3. Click on the "download" button

Here are my six tasks and questions for this week. As usual, a link to the Jupyter notebook I used to solve these is at the bottom of the post.

Download the CSV data. Read it into a data frame using both the PyArrow engine and for all of the backend dtypes (i.e., instead of NumPy).

Before doing anything else, let’s load Pandas:

import pandas as pd

I asked you to load the CSV file into a data frame, which we can do with “read_csv”:

filename = '2021_Final_Assisted_Reproductive_Technology__ART__Summary_20240228.csv'

df = pd.read_csv(filename)

This works just fine, but it takes a long time. That’s because by default, Pandas uses its own engine for reading CSV files, and uses NumPy for back-end storage. Both have served the Pandas community for a long time, but they are slowly being phased out in favor of PyArrow (https://pypi.org/project/pyarrow/), Python bindings for the Apache Arrow open-source project (https://arrow.apache.org/). Arrow is an in-memory database optimized for the sorts of things that Pandas does, which means that it is often faster than the default implementations. It also takes less memory, partly thanks to its native use of strings. Pandas is already warning us that version 3.0 will require PyArrow.

Today, we can use PyArrow in two different ways:

  1. We can specify that Pandas use it for parsing our CSV file, and
  2. We can use it instead of NumPy on the back end.

We can do both by passing the “engine” and “dtype_backend” keyword arguments to “read_csv”:

filename = '2021_Final_Assisted_Reproductive_Technology__ART__Summary_20240228.csv'

df = pd.read_csv(filename, engine='pyarrow', dtype_backend='pyarrow')

I compared the time it took to with the default settings and with the PyArrow engine and dtypes. On average, the default settings took 201 ms. By contrast, using PyArrow took only 12.9 ms on average.

Moreover, using the default NumPy dtypes consumed 88,423,682 bytes. Using PyArrow, we only needed 29,528,256 bytes, or less than half.

In other words: PyArrow loaded our data by a factor of nearly 20x the speed, and used half the memory, without any additional configuration or customization. PyArrow isn’t always faster than NumPy, and it doesn’t handle all of the functionality that NumPy does. But I’d argue that it’s worth using it wherever you can.

How many IVF clinics were surveyed for this data? How many questions were asked of each clinic?

It took me some time to understand how this data set was structured. I originally thought that each row of the CSV file represented a single fertility clinic, and that the columns were their answers to various questions. That’s not exactly the case; each question was listed in a separate row, which means (in theory) that if there were m clinics and n questions, then the data set has m*n rows.

The question text is in the “Question” column, while the answer is actually located in two columns, one called “Data_Value” which contains whatever the response was in string dtypes, and one called “data_value_num” which contains numeric values only as float dtypes. Because we’re using PyArrow, string values are “string[pyarrow]” and float values are “double[pyarrow]”.

Finding the number of clinics surveyed thus required that we get the unique values associated with the clinics.

Rows also contained “Breakout groups” and “Breakout” columns, which you can think of a form of multi-index inside of the CSV file. The breakout groups were the age of the patient, the egg/embryo type, and a generic yes/no group. The values themselves were age ranges, and types of eggs and embryos, among others. Because of the breakout groups and values, we have many more rows than the m*n calculation I mentioned above; again, if we think of the breakout groups as a multi-index, we can understand how the clinics were answering the survey questions not for all of their patients, but for specific groups, to give the CDC more precise data.

With a bit more work, we could have broken down the data in all sorts of ways — including the effectiveness of IVF for different age ranges and the effectiveness of different techniques. However, I decided that it was complex enough to deal with the questions, and focused on those for my questions this week.

And thus, how can we find out how many different clinics were surveyed? I decided to look at the “FacilityName” column, and to ask Pandas for the number of unique values it contained, using the “nunique” method:

df['FacilityName'].nunique()

I got a result of 453. Another way to do this would be to invoke the “unique” method:

df['FacilityName'].unique()

This returns a NumPy array of the unique values. There’s a third method we can try, “drop_duplicates”, which returns a Pandas series of the values:

df['FacilityName'].drop_duplicates()

Which of these does the job fastest? I again ran “%timeit” in Jupyter, and found:

We again see that while Pandas offers many ways to accomplish the same (or similar) goal, choosing well can have a big impact on your query performance.

As for how many questions were asked, I ran the following:

df['Question'].nunique()

We see that the survey asked clinics to answer 45 questions.

Which state had the largest number of clinics in this study, and which had the smallest number?

After reading about the Alabama court ruling, I was curious about how many clinics in the US even offer IVF. I thus asked you to find out how many clinics are in each state, and then to show the smallest and largest numbers.

First, I kept only the two columns that I need in order to make this determination, “LocationAbbr” (the two-letter state abbreviation) and “FacilityName” (the clinic’s name):

(
    df
    [['LocationAbbr', 'FacilityName']]
)

I’ll now use “drop_duplicates” to get unique rows back from our data frame:

(
    df
    [['LocationAbbr', 'FacilityName']]
    .drop_duplicates()
)

Why didn’t I use either “unique” or “nunique”? Because they only work on a series, and I have a data frame here.

Next, I ask for only the “LocationAbbr” column:

(
    df
    [['LocationAbbr', 'FacilityName']]
    .drop_duplicates()
    [['LocationAbbr']]
)

Now that I have the locations, I can count them with “value_counts”:

(
    df
    [['LocationAbbr', 'FacilityName']]
    .drop_duplicates()
    [['LocationAbbr']]
    .value_counts()
)

That tells me how often each state appears. I could call “head” to get the most common (since they’re at the top of the result) or “tail” to get the least common (since they’re at the bottom), but I can also use “agg” to get both:

(
    df
    [['LocationAbbr', 'FacilityName']]
    .drop_duplicates()
    [['LocationAbbr']]
    .value_counts()
    .agg(['head', 'tail'])
)

The result:

              head  tail
LocationAbbr            
CA            78.0   NaN
NY            45.0   NaN
TX            43.0   NaN
FL            26.0   NaN
IL            24.0   NaN
MT             NaN   1.0
ME             NaN   1.0
ID             NaN   1.0
AR             NaN   1.0
AK             NaN   1.0

Notice how we get NaN values for wherever there are no values. Also notice that because NaN is a float value, it forced the two columns to be floats.

Also notice that yes, the states in the “tail” output all have values of 1.0. They aren’t the only states with 1 clinic in the survey — but that’s all “tail” will show us.

Not surprisingly, the largest states have the greatest number of clinics, with California way out in front with 78, followed by New York, Texas, and Florida.

Create a two-column data frame in which the index is the clinic name, and the columns contain the mean (a) number of retrievals and (b) number of transfers done by that clinic.

Wrestling a data frame into shape can be difficult — but since restructuring our data often makes it easier to work with, it’s also a common task. Here, I asked you to take the entire data set, grab the answers to two of the questions, and then rejigger the data frame so that we can make some sense of the data.

First, we want to see only those rows from the “Number of retrievals” and “Number of transfers” questions. To get these, I used “loc” along with a lambda expression that returns True/False for each row in the data frame, in which I use “isin” to tell me whether the question is one of the two I’m looking for:

(
    df
    .loc[lambda df_: df_['Question'].isin(['Number of retrievals', 
            'Number of transfers'])]
)

Now that I have pared my data frame down to rows for these two questions, I pare it down further, only looking at three columns — the facility name, the question, and the numeric answer:

(
    df
    .loc[lambda df_: df_['Question'].isin(['Number of retrievals', 
            'Number of transfers'])]
    [['FacilityName', 'Question', 'data_value_num']]
)

I then run “groupby”, grouping on two columns (“FacilityName” and “Question”). For each unique combination of values in these columns, I calculate the mean of “data_value_num”:

(
    df
    .loc[lambda df_: df_['Question'].isin(['Number of retrievals', 
            'Number of transfers'])]
    [['FacilityName', 'Question', 'data_value_num']]
    .groupby(['FacilityName', 'Question'])
    ['data_value_num'].mean()
)

The good news is that this gives me a series back, one with a 2-level multi-index whose values show the mean number of retrievals or transfers.

The bad news is that this isn’t what I asked for: I wanted a data frame with the facilities in the index and the questions in the columns. That’s where “unstack” comes into play; it moves part of a multi-index into column names. I just had to specify that I want to take the inner part of the multi-index (i.e., level=1), and it worked:

(
    df
    .loc[lambda df_: df_['Question'].isin(['Number of retrievals', 
            'Number of transfers'])]
    [['FacilityName', 'Question', 'data_value_num']]
    .groupby(['FacilityName', 'Question'])
    ['data_value_num'].mean()
    .unstack(level=1)
)

The result is a data frame with 453 rows (one for each clinic) and two columns (one for each question).

How many total retrievals and transfers were done in Alabama?

To answer this, we’ll need to find rows where:

  1. The question is either “Number of retrievals” or “Number of transfers” and
  2. The LocationAbbr column is “AL”

The first part is identical to what we did in the previous question. But how can we use the boolean series that we got back from that, and combine it with the requirement that the location be in Alabama?

We can use the “&” operator, which takes two boolean series with the same index, and returns a new boolean series with that same index — one in which each value is False, unless both of the original series are True. Yes, it’s the Pandas “and” operator. We can’t use the standard Python “and” operator, because it works with individual boolean values, and here we have series (aka vectors) of booleans. Pandas, like NumPy, took advantage of operator overloading and the “&” bitwise operator.

I’m also only interested in the “Question” and “data_value_num” columns, so I use the two-argument version of “loc”, selecting columns in that second argument:

(
    df
    .loc[(df['LocationAbbr'] == 'AL') & (df['Question'].isin(['Number of retrievals', 'Number of transfers'])), 
    ['Question','data_value_num']]
)

With only Alabama rows in the data frame, and only the two questions in our “Questions” column, I can now calculate how many total retrievals and transfers were done in Alabama in 2021:

(
    df
    .loc[(df['LocationAbbr'] == 'AL') & (df['Question'].isin(['Number of retrievals', 'Number of transfers'])), 
    ['Question','data_value_num']]
    .groupby('Question').sum()
)

I got results of 486 retrievals and 521 transfers.

Show, for each state, the mean percentage of live-birth deliveries after 1, 1 or 2, or all intended retrievals. The index should be the states, and the columns should contains the mean answers to the questions that start with "Percentage of new patients having live-birth deliveries after".

Finally, I wanted to see a data frame in which:

This sounds suspiciously like a pivot table, no? Remember that a pivot table is basically a two-dimensional groupby: We take two categorical columns and one numeric column, and then run “mean” on each two-category intersection.

To do this, we’ll first select the three columns we’ll need:

(
    df
    [['LocationAbbr', 'Question', 'data_value_num']]    
)

To get only those rows with the questions we can use “str.contains”, passing the regex=True option to use a regular expression:

(
    df
    [['LocationAbbr', 'Question', 'data_value_num']]
    .loc[lambda df_: df_['Question'].str.contains('^Percentage of new patients having live-birth deliveries after', regex=True)]       
)

With this in place, we can call “pivot_table”:

(
    df
    [['LocationAbbr', 'Question', 'data_value_num']]
    .loc[lambda df_: df_['Question'].str.contains('^Percentage of new patients having live-birth deliveries after', regex=True)]
    .pivot_table(index='LocationAbbr', columns='Question', values='data_value_num')
    
)

That gives us the result we want, a data frame with 51 rows (for each state + DC) and three columns. I don’t know about you, but this data tells me that IVF is rather effective; a very high proportion of patients ended up having live births. Which is why it’s such a widely used technique, and why the Alabama court’s decision set off such shock waves.

What do you think? What other analysis could/should we have done? Leave a comment!

Meanwhile, here is my Jupyter notebook: https://drive.google.com/file/d/1hPmwc_L_xUVjgfb5va00iLNWz7e-TyDh/view?usp=sharing

I’ll be back next week with more questions and puzzles in data analytics with Pandas.

Reuven