Module 6: Reproducible pipeline

PUBHLT 0411 · Python for Public Health Data Analysis

Soumik Purkayastha

Dept. of Biostatistics and Health Data Science, University of Pittsburgh

November 17, 2026

A finding that contradicts what everyone knows.

Firearm death is understood as an urban problem, and Pennsylvania’s data disagrees

Firearm injury is commonly described as a problem of cities. Pennsylvania’s county-level data, correctly assembled, reports a higher firearm death rate outside metropolitan counties than inside them.

The question this session opens

Which Pennsylvania counties carry the highest burden of firearm death — and does the answer depend on how the analysis handles the counties that report nothing?

This session establishes the pipeline as a single reproducible object: import, clean, join, group, summarise, visualise, report — one document another person can run.

Required from Modules 4 and 5: pd.read_csv(), .dropna(), .merge(), .groupby(), .agg(), and plt.subplots().

Ask the room first, before any code is on screen. Where in Pennsylvania would you expect the firearm death rate to be highest? Take three answers. Somebody says Philadelphia, somebody says Pittsburgh, and that is the honest expectation — it is what the news covers and it is not a stupid answer.

Then tell them the data disagrees, and do not tell them by how much yet. Today they build the analysis that shows it.

Here is what makes this session different from every other one. In five modules the finding confirmed what a public-health student already suspected. Poverty tracks infant mortality. Smoking tracks asthma. Those were real, and they were reassuring. This one is not reassuring, and it is the only finding in the course designed to surprise them.

And the reason it surprises is the part worth teaching. It is not that everyone was wrong about cities. It is that the obvious way to compute this number gives the opposite answer, and the obvious way is wrong for a reason they can see and fix. That is the whole session.

One thing to say out loud early, because it sets the frame for the next seventy-five minutes. This is injury prevention. Roughly six in ten of these deaths are suicides, which is a mental health and means-access story, not a crime story. This course treats it that way, and the writeup they turn in should too.

This session is one working session, and the project starts in it

First twenty minutes	The pipeline as one object, and the one error that reverses this finding
The rest of the session	Lab — build the analysis, end to end, on the real data
What you leave with	A document that runs top to bottom, and the paragraph that reports it

The pipeline is seven stages, and the whole session is one pass through it

Stage	Verb	Module
Import	pd.read_csv()	4
Clean	.dropna(), .fillna(), .astype()	4
Join	.merge()	4
Group	.groupby()	4
Summarise	.agg(), derived columns	4
Visualise	plt.subplots(), ax.bar()	5
Report	prose, with the numbers in it	6

Every stage has already been taught. The new skill is running all seven as one object that another person can execute and get the same answer.

Point at the right-hand column, because it is the reassuring part. There is exactly one new row in this table, and it is the last one. Everything above it they already did in the last two modules.

What is new is the seam between the stages. Individually they can each merge and each group. The thing they have never done is hold all seven in one document in an order that survives being run by somebody else, on a different machine, three weeks later. That is a genuinely different skill and it is the one employers actually notice.

The word to keep saying today is pipeline, singular. Not “my analysis code.” One object, with an input at one end and a claim at the other.

The data is two files, and the join is the reason there are two

pa_firearm_county_year.csv has 402 rows — 67 counties × 6 years, 2019–2024. It carries deaths and it does not carry metro status. pa_counties.csv carries metro status and it does not carry deaths.

Source

CDC, Mapping Injury, Overdose, and Violence — County (psx4-wq38); Census CO-EST2023 and the NCHS Urban–Rural Classification, 2013.

This is the honest motivation for the join, and it is the same one from Module 4. The question is about metro against nonmetro. The word “metro” appears in exactly one of these two files, and the deaths appear in the other. There is no way to answer the question without putting them together, which is why merge is not a party trick.

Worth naming the classification, since a student will ask what counts as metro. It is the NCHS urban–rural scheme, six levels collapsed to two, and it is assigned by the federal government rather than by the instructor. Thirty-seven metro counties, thirty nonmetro.

Forty percent of the cells report no number, and the reason is protective

CDC suppresses a county-year cell when the count is small enough that publishing it could identify a decedent. 159 of 402 cells carry no number. In 2023, 25 of 67 counties report nothing.

Suppression is not random, and that is the whole problem

Cells are suppressed because the count is small. Small counts occur in small counties, and small counties are disproportionately nonmetro — 19 of 30 nonmetro counties are suppressed in 2023, against 6 of 37 metro counties.

Take a moment on why suppression exists, because students arrive assuming missing data is sloppiness. It is the opposite of sloppiness. If a county of four thousand people had one firearm death last year and the agency published it, somebody reading the local paper could work out who. The blank is a deliberate act of protection, and the right response to it is not irritation.

Now the sentence that matters for today. Missingness here is not a coin flip. It is caused by the very thing being measured — small counts — and small counts live in small rural counties. Nineteen of thirty nonmetro counties versus six of thirty-seven metro. That is a structural pattern, not noise.

Hold that thought for two slides, because it is about to eat the analysis alive.

The pooled rate is total deaths over total population

A pooled rate sums the numerators and sums the denominators, then divides once. It is not the average of the county rates.

Walk the four lines slowly, because this is the shape of every pipeline they will ever write. Filter to a year. Drop the rows with no number. Join on the grouping variable. Group and aggregate. Four verbs, all from Module 4.

The one line to defend is the dropna, and it is on the next slide, so resist explaining it here. Let them see the output first.

The distinction at the bottom is worth thirty seconds. Pooling is not averaging. If the sixty-seven county rates are averaged with equal weight, Cameron County with four thousand people counts as much as Philadelphia with a million and a half. Pooling weights by population, which is what a state-level rate means. Both are legitimate and they answer different questions; today is pooled.

The same analysis without the .dropna() reverses the finding

2023, pooled	Correct	Without .dropna()
Metro	14.1	13.6
Nonmetro	18.5	8.2

.sum() skips missing values, and the population column does not. A suppressed county therefore contributes its full population to the denominator and zero deaths to the numerator, which drags the group rate downward — hardest where suppression concentrates.

This error raises no exception

Both versions run. Both produce a clean table. One says nonmetro rates are lower.

This is the most important slide in the module. Slow all the way down.

Ask them to read the second column out loud. Eight point two against thirteen point six. Nonmetro looks safer, by a lot, and it confirms exactly what the room guessed twenty minutes ago. If a student handed that table in, it would look completely reasonable. The code is clean. Nothing errored. The numbers have the right number of decimal places.

And it is wrong, for a reason that fits in one sentence: sum skips NaN and the denominator does not. Nineteen nonmetro counties put their entire population into the bottom of that fraction and contributed no deaths to the top, because their deaths were suppressed for privacy. The rate is not measuring firearm death. It is measuring how much nonmetro population got parked in a denominator with nothing above it.

The fix is one line and they have known it since Module 4. Drop the rows with no number before the join, so a county that reports nothing is out of both halves of the fraction.

Then the sentence that is really about the rest of their careers. This error is silent. There is no exception, no warning, no red text. In Module 5 the lesson was that a figure can run and lie. Today it is that a table can run and lie, and the only defence is knowing what your denominator contains. Say it exactly like that.

The suicide share explains the gap

Among counties reporting a death count, 62.1% of nonmetro firearm deaths are suicides, against 54.7% in metro counties.

The definition is part of the finding

Suppressed suicide counts are treated as zero here, across the counties that report a death total. Restricting instead to counties reporting both columns gives 83.7% for nonmetro — a different, much smaller denominator. The number moves when the definition moves, so the definition is stated with the number.

This is the answer to the “but why” that the previous slide raises. A higher nonmetro rate is a finding. It is not yet an explanation, and a finding without a mechanism is where a lot of student writeups stop.

Sixty-two percent against fifty-five. The nonmetro burden is more heavily suicide, and that is consistent with everything known about firearm injury in rural America — higher household firearm ownership, greater physical isolation, thinner mental health infrastructure, and longer distances to emergency care. Suicide attempts with a firearm are lethal in roughly eighty-five percent of cases, against a few percent for most other means, so access to a firearm converts a crisis into a death at a rate nothing else matches.

Now the callout, which is a methods lesson wearing a public health hat. There are two defensible ways to compute that share and they differ by twenty-two points. Not because anyone cheated — because the suicide column is suppressed more often than the death column, so restricting to rows with both quietly throws away counties. Whichever one is chosen, it gets stated in the writeup. An unlabelled 62% is not a number, it is a rumour.

The finding holds across six years, so it is not one unusual year

Figure 1. Line chart of pooled firearm death rate per 100,000 for reporting Pennsylvania counties, 2019 to 2024, with two series distinguished by both colour and marker shape. Nonmetro counties, drawn in dark gold with triangular markers, run at 18.0 in 2019, dip to 13.4 in 2020, then rise to 16.3, 16.1, and 18.5 before falling to 15.7 in 2024. Metro counties, drawn in royal blue with circular markers, rise from 11.6 in 2019 to a peak of 15.3 in 2022, then fall to 14.1 and 12.4. The nonmetro line sits above the metro line in five of the six years; the exception is 2020.

This is the reproducibility argument made visually. One year could be a fluke, and a careful reader is entitled to ask whether 2023 was cherry-picked. Five of six years say the same thing, and the exception is 2020, which is the pandemic year that behaves strangely in every dataset in this course.

Point at the marker shapes deliberately. Triangles and circles, not just gold and blue. From Module 5: those two colours are each perfectly legible against white and they do not reliably separate from each other — two point oh two to one. Colour is never the only channel, and that rule does not relax because the figure is going in a report.

The zero baseline matters here too. Setting the y-limit to start at zero rather than at thirteen is what keeps a four-point gap looking like a four-point gap.

Mention the counties-reporting-per-year caveat if anyone is sharp enough to ask: the nonmetro line is built from seven to eleven counties depending on the year, so it moves more than the metro line. That is real instability, and it is the honest thing to put in the writeup rather than hide.

Notebook hygiene is what makes the result reproducible

Rule	Reason
Imports in the first cell	A notebook that only runs after scrolling up is broken
Restart and run all, before submitting	Out-of-order state is invisible until it is fatal
Relative paths — data/file.csv	/Users/soumik/Desktop/... runs on exactly one machine
np.random.seed(...) if anything is random	Same input, same output, every time
No manual edits to output	A number typed by hand is not a result

The failure this prevents

A notebook where cell 12 was run, then cell 4 was edited, holds a value in memory that no longer matches any code on screen. The document is the record, not the memory.

This is the section that sounds like housekeeping and is actually the difference between an analysis and a rumour.

Tell them the restart-and-run-all story, because everyone has a version of it. A student builds something over four hours, it works, they close the laptop. They reopen it the next morning and nothing runs, because half the variables they depended on were defined in cells they later deleted or edited. The notebook never told them, because the values were sitting in memory the whole time. Kernel restart, run all, from the top — before submitting anything. Every time. It takes thirty seconds and it is the single highest-value habit in this course.

The absolute-path one is worth a live demonstration if there is time. A path starting with slash Users slash somebody’s name runs on precisely one computer on Earth. Their project gets opened by a grader on a different machine, and a data-slash-filename path just works.

The last row is the ethics row, so do not skip it. If the code prints 18.4 and the writeup says 18.5 because that looked tidier, that is fabrication. Small, harmless-feeling, and still fabrication. The number in the prose comes from the cell above it.

Uncertainty is reported in words, with the counts that produced it

The nonmetro rate is built from 11 counties and 116 deaths. The metro rate is built from 31 counties and 1,567 deaths. A rate computed from a small number of events moves a great deal when one event is added or removed.

Written honestly	Written badly
“18.5 per 100,000, from 116 deaths across 11 reporting counties”	“18.5 per 100,000”
“Suppression removes 19 of 30 nonmetro counties”	“Nonmetro counties: 18.5”
“The gap holds in five of six years”	“Nonmetro rates are higher”

This course does not compute p-values

Uncertainty is described: how many events, how many counties, how much the estimate moves between years, and what suppression removed. A p-value is not required to report an estimate honestly, and it is not a substitute for reporting the counts.

This is the §8 slide and it is the one an outside reviewer would flag, so say the reasoning out loud rather than hiding it.

They have not had inference. No t-tests, no confidence intervals, no p-values. That is a deliberate choice about what fits in six weeks, and here is the defence: a p-value they cannot interpret is worth less than a count they can. Telling a reader the estimate rests on 116 deaths in 11 counties conveys the fragility directly, in a form nobody misreads. A star next to a number does not.

Juniata County is the concrete example and it is a good one. Ten deaths, twenty-three thousand people, forty-three per hundred thousand — the highest county rate in the state that year. Ask the room what happens to that number if two of those ten deaths had happened in December of the previous year. It drops to thirty-four. Same county, same reality, different bin. That is what an unstable rate means, and no student needs a distributional assumption to feel it.

For their projects the rule is simple: never report a rate without the count it came from.

The report is the deliverable, and it has a fixed shape

Move	One or two sentences
The question	What was asked, of which population, in which years
The data	Source, coverage, and what is missing and why
The method	The pipeline in prose — filtered, dropped, joined, pooled
The finding	The numbers, with their denominators
The limit	What this design cannot establish

Integrating computer output

Output is evidence inside a sentence, never a pasted block left to speak for itself. ✗ A bare table dropped into the document. ✓ “Among the 42 counties reporting a 2023 count, the pooled nonmetro rate was 18.5 per 100,000 against 14.1 in metro counties (Table 1).”

Give them the shape explicitly, because “write it up” is useless advice and this table is not. Five moves. A paragraph each, or less. That is a report.

The move students skip is the last one, and it is the one that separates a public-health writeup from a blog post. What this design cannot establish: these are county aggregates, so a pattern across counties describes counties and not the individuals in them — that is the ecological fallacy, and it applies here exactly as it did to smoking and asthma last week. Metro status is not a cause of anything. It is a label that travels with firearm ownership, isolation, and access to care.

On the callout: the commonest failure in the drafts is a screenshot of a DataFrame sitting alone in the middle of a page, doing no work. A table is an exhibit. A sentence points at the exhibit and says what to see in it. Show them the good version and read it aloud — it names the coverage, gives both numbers, and cites the table, in one breath.

One more, and be firm because it comes up every year: generative AI is not permitted in this course. The writeup is assessed as their reasoning about their own analysis.

The finding, stated as it would appear in a report

Among the 42 Pennsylvania counties reporting a 2023 firearm death count, the pooled death rate was 18.5 per 100,000 in nonmetropolitan counties and 14.1 per 100,000 in metropolitan counties. Suicides account for 62.1% of nonmetropolitan firearm deaths and 54.7% of metropolitan ones. Because CDC suppresses small counts, 25 of 67 counties report no 2023 total, and suppression is concentrated in nonmetropolitan counties (19 of 30). Counties reporting no total are excluded from both the numerator and the denominator; retaining their population while summing only reported deaths reverses the comparison to 8.2 against 13.6. The nonmetropolitan rate exceeds the metropolitan rate in five of the six years from 2019 to 2024.

What this paragraph does

It states the numbers, the coverage, the definition, the sensitivity, and the consistency check. It makes no causal claim.

Read this paragraph out loud, slowly, and then say: that is the whole assignment. Not the notebook — this. Six sentences, every one carrying a number, no adjectives doing any work.

Have them count what is in it. The estimate. The denominator. The definitional choice. The sensitivity analysis, which is the sentence about 8.2 and 13.6 — that is the alternative specification, named and reported rather than buried. And the robustness, which is five of six years. That is a professional paragraph, and no student in the room could not have written it by the end of today.

Then point at what is not in it. No “shockingly.” No “this proves.” No mechanism asserted without evidence. The suicide share is offered as consistent with an explanation, and the paragraph stops there, because that is where the data stops.

This is the target for the lab. Their project paragraph will have a different question in it, and the same five moves.

Lab

The remainder of the session is build time. Nine exercises assemble the pipeline in the order it runs, and the last one produces the paragraph.

1–4	Import, inspect, clean, and join — including the suppression decision
5	Question only, no scaffold — group and summarise, for the headline finding
6	The silent error, reproduced deliberately and repaired
7–8	Question only — the figure and its description
9	Question only — write the finding as a paragraph

Four of these supply no code

Exercises 5, 7, 8, and 9 state a question and the deliverable, and the cell is empty. The project will not supply a scaffold either.

This is the project skeleton

The final project replaces the question and the dataset. The seven stages do not change.

Set them going and then circulate rather than lecturing over the top of it.

Warn them about the shape before they start. The first four exercises hand them a scaffold with gaps. Exercise 5 does not — it asks a question in English and gives them an empty cell, and so do 7, 8, and 9. That jump is deliberate and it is worth naming, because a student who hits an empty cell unprepared assumes something is broken. Tell them the hint is there and that using it costs nothing.

Exercise 5 is where most hands will go up, since it is the first blank page. The move that unsticks people is asking what the table should look like when it is finished — how many rows, what the columns are called. Once they can describe the output they can usually write the groupby.

Exercise 6 is the other one to watch. It hands them the wrong pipeline, running cleanly, and asks what is wrong with it — and a good number of students will stare at correct-looking code for several minutes before the denominator occurs to them. That struggle is the point, so do not rescue anyone too early. Ask them what is in the bottom of the fraction and walk away.

For anyone who finishes early, the extension is in the final callout: run the same pipeline grouped by region instead of metro status, six groups instead of two, and see whether the story survives a finer cut. It mostly does, and the counties get thin enough that the instability lesson lands a second time on its own.

Remind them at about the twenty-minute mark to restart and run all. Some of them will already have out-of-order state, and finding that now is much cheaper than finding it in December.

Exercise 1 — Import and inspect

Exercise 1 — Import and inspect

Worked example. A pipeline opens by loading its inputs and confirming their shape. An analysis that never checks what arrived is an analysis that trusts a file it has not read.

Exercise 1 — Import and inspect — your turn

Your turn. Load data/pa_counties.csv into county_reference. Print its shape, how many counties fall into each metro category, and confirm the file has one row per county.

# Four blanks: read_csv, shape, value_counts, nunique.
#
# The last line is the one worth thinking about. It asks whether the
# number of distinct county names equals the number of rows — if a county
# appeared twice, every rate you compute later would be wrong, and nothing
# would raise an error. Checking it at import costs one line.

import pandas as pd county_reference = pd.read_csv("data/pa_counties.csv") print(county_reference.shape) print(county_reference["metro"].value_counts()) print("one row per county:", county_reference["county"].nunique() == len(county_reference))

import pandas as pd

county_reference = pd.read_csv("data/pa_counties.csv")

print(county_reference.shape)
print(county_reference["metro"].value_counts())
print("one row per county:", county_reference["county"].nunique() == len(county_reference))

Sixty-seven, every time

Pennsylvania has 67 counties. A county-level table that does not have 67 rows has either lost some or gained duplicates, and finding that at import costs nothing.

Exercise 2 — Filter to one year

Exercise 2 — Filter to one year

Worked example. The file holds six years stacked in one column. A question about 2023 requires selecting those rows before anything else happens.

Exercise 2 — Filter to one year — your turn

Your turn. Build f2023, the 2023 rows only. Print its shape, how many firearm_deaths values are missing, and what fraction of counties that leaves reporting.

# Four blanks: firearm, firearm, ==, isna.
#
# A single = would try to assign, which is not what a filter does. You are
# asking each row a yes/no question, exactly as in Module 1.
#
# Note what the last line reports. Not a count — a proportion. "25 missing"
# means nothing without "of 67", and that fraction is the first thing a
# reader of your report will want.

f2023 = firearm[firearm["year"] == 2023] n_missing = f2023["firearm_deaths"].isna().sum() n_reporting = len(f2023) - n_missing print(f2023.shape) print(n_missing, "suppressed in 2023") print("reporting:", n_reporting, "of", len(f2023), f"({n_reporting / len(f2023):.0%})")

f2023 = firearm[firearm["year"] == 2023]

n_missing = f2023["firearm_deaths"].isna().sum()
n_reporting = len(f2023) - n_missing

print(f2023.shape)
print(n_missing, "suppressed in 2023")
print("reporting:", n_reporting, "of", len(f2023),
      f"({n_reporting / len(f2023):.0%})")

Exercise 3 — Clean, and decide about the blanks

Exercise 3 — Clean, and decide about the blanks

Worked example. .dropna(subset=[...]) removes rows missing a value in the named column and leaves every other column alone. Naming the column matters: dropping on the whole frame would also discard counties whose homicide count is suppressed but whose total is known.

Exercise 3 — Clean, and decide about the blanks — your turn

Your turn. Build reporting: the 2023 counties that report a firearm_deaths total, dropping on that column only. Then show what the careless version costs, by naming the counties it would discard even though their death total is known.

# Four blanks: dropna, subset, dropna, sum.
#
# The `lost` line is the point of the exercise. It subtracts the careless
# result from the careful one, so what remains is the counties whose death
# total is perfectly good but whose homicide or suicide cell happens to be
# suppressed. A bare .dropna() silently discards all 23 of them.
#
# Read the names it prints. Bucks, Chester, Westmoreland, Dauphin — these
# are among the largest counties in the state. Losing them is not a small
# rural rounding error; it would gut the metro group.

reporting = f2023.dropna(subset=["firearm_deaths"]) lost = set(reporting["county"]) - set(f2023.dropna()["county"]) print(len(reporting), "counties report a 2023 total") print(reporting["firearm_deaths"].sum(), "deaths among them") print(len(lost), "counties a bare dropna() would have thrown away:") print(sorted(lost))

reporting = f2023.dropna(subset=["firearm_deaths"])

lost = set(reporting["county"]) - set(f2023.dropna()["county"])

print(len(reporting), "counties report a 2023 total")
print(reporting["firearm_deaths"].sum(), "deaths among them")
print(len(lost), "counties a bare dropna() would have thrown away:")
print(sorted(lost))

Twenty-five counties just left the analysis

They left for a defensible reason, and the report must say so. A cleaning step that is not described in the writeup is a cleaning step the reader cannot check.

Exercise 4 — Join

Exercise 4 — Join

Worked example. .merge() attaches columns from a second table by matching on a shared key. how="left" keeps every row of the left table whether or not it finds a match, which makes an unmatched key visible instead of silent.

Exercise 4 — Join — your turn

Your turn. Merge reporting with the county, metro, and population_2023 columns of county_reference, matching on county. Confirm the row count is unchanged and that nothing failed to match.

# Five blanks: merge, on, left, ==, isna.
#
# The row-count check catches a failure the missing-value check cannot. If
# the right-hand table had two rows for some county, the merge would
# silently DUPLICATE the left-hand row and that county's deaths would be
# counted twice. A join can lose rows and it can invent them; check both.
#
# `on="county"` says the two tables are lined up by county name. If these
# files spelled a county differently — the real "Mc Kean" / "McKean" bug
# from Module 4 — that county would arrive with a missing metro value,
# which is exactly what the second print statement is checking for.

before = len(reporting) joined = reporting.merge( county_reference[["county", "metro", "population_2023"]], on="county", how="left" ) print(joined.shape) print("rows preserved:", len(joined) == before) print(joined["metro"].isna().sum(), "rows failed to match")

before = len(reporting)

joined = reporting.merge(
    county_reference[["county", "metro", "population_2023"]],
    on="county",
    how="left"
)

print(joined.shape)
print("rows preserved:", len(joined) == before)
print(joined["metro"].isna().sum(), "rows failed to match")

Exercise 5 — Group and summarise

Exercise 5 — Group and summarise

Worked example. .groupby().agg() collapses many rows into one row per group, naming each output column and the operation that fills it. "size" counts rows; "sum" totals a column.

Exercise 5 — Group and summarise — your turn

Your turn — no scaffold. Answer this question from joined:

Do Pennsylvanians living outside metropolitan counties die by firearm at a higher rate than those living in them?

Produce a table called summary, one row per metro group, reporting for each group: how many counties it contains, its total firearm deaths, its total population, and its firearm death rate per 100,000 people. Print it.

# Two verbs and one arithmetic line, all of which you have used before.
#
# Group `joined` by the metro column, and inside .agg() name each output
# column and the (column, operation) pair that fills it. "size" counts the
# rows in a group; "sum" totals a column. That gets you counties, deaths,
# and population.
#
# Then the rate. A rate is events divided by the population at risk, times
# 100,000 — which is what makes it "per 100,000" rather than a proportion.
# Add it as a new column on the grouped table, not on `joined`.
#
# Name the columns exactly: counties, deaths, population, rate_per_100k.

summary = joined.groupby("metro").agg( counties=("county", "size"), deaths=("firearm_deaths", "sum"), population=("population_2023", "sum"), ) summary["rate_per_100k"] = summary["deaths"] / summary["population"] * 100000 print(summary.round(1))

summary = joined.groupby("metro").agg(
    counties=("county", "size"),
    deaths=("firearm_deaths", "sum"),
    population=("population_2023", "sum"),
)

summary["rate_per_100k"] = summary["deaths"] / summary["population"] * 100000
print(summary.round(1))

That is the finding

18.5 against 14.1. The rate outside metropolitan Pennsylvania is higher, and it was produced by five ordinary lines of pandas.

Exercise 6 — The error that raises nothing

Exercise 6 — The error that raises nothing

Worked example. The pipeline below is identical to Exercise 5 except for one omission. It runs without complaint and returns the opposite finding.

The suppressed counties were never removed. .sum() skips their missing deaths; the population column does not skip their people.

Exercise 6 — The error that raises nothing — your turn

Your turn. Repair it, then measure the damage the broken version did: how many people were sitting in the nonmetro denominator with no deaths above them, and how far that pushed the nonmetro rate.

# Five blanks: dropna, subset, population, -, -.
#
# Think about what belongs in a denominator. A rate's denominator is the
# population a counted event could have happened to. If a county's deaths
# were never published, its people cannot sit in the bottom of the fraction
# while nothing of theirs is in the top.
#
# The last two lines quantify the error rather than just correcting it.
# That number — how far wrong the careless version was — is what belongs in
# a writeup's sensitivity sentence.

repaired = f2023.dropna(subset=["firearm_deaths"]).merge( county_reference[["county", "metro", "population_2023"]], on="county", how="left" ) fixed = repaired.groupby("metro").agg( deaths=("firearm_deaths", "sum"), population=("population_2023", "sum"), ) fixed["rate_per_100k"] = fixed["deaths"] / fixed["population"] * 100000 phantom = (broken_summary.loc["Nonmetro", "population"] - fixed.loc["Nonmetro", "population"]) understated_by = (fixed.loc["Nonmetro", "rate_per_100k"] - broken_summary.loc["Nonmetro", "rate_per_100k"]) print(fixed.round(1)) print(f"{phantom:,} nonmetro residents were in the denominator with no deaths above them") print(f"the broken rate understated nonmetro by {understated_by:.1f} per 100,000")

repaired = f2023.dropna(subset=["firearm_deaths"]).merge(
    county_reference[["county", "metro", "population_2023"]], on="county", how="left"
)

fixed = repaired.groupby("metro").agg(
    deaths=("firearm_deaths", "sum"),
    population=("population_2023", "sum"),
)
fixed["rate_per_100k"] = fixed["deaths"] / fixed["population"] * 100000

phantom = (broken_summary.loc["Nonmetro", "population"]
           - fixed.loc["Nonmetro", "population"])
understated_by = (fixed.loc["Nonmetro", "rate_per_100k"]
                  - broken_summary.loc["Nonmetro", "rate_per_100k"])

print(fixed.round(1))
print(f"{phantom:,} nonmetro residents were in the denominator with no deaths above them")
print(f"the broken rate understated nonmetro by {understated_by:.1f} per 100,000")

Running is not checking

The broken version produced a clean table with plausible numbers and no warning. Every result in a report needs a reason to be believed that is stronger than the absence of an error message.

Exercise 7 — Visualise

Exercise 7 — Visualise

Worked example. A two-bar comparison is the right figure for a two-group finding. The zero baseline is what keeps the comparison honest, and the value labels remove the need to read heights against a gridline.

Figure 2. Bar chart of the pooled 2023 firearm death rate for the 42 reporting Pennsylvania counties. The metro bar, in royal blue, reaches 14.1 deaths per 100,000. The nonmetro bar, in dark gold, reaches 18.5. Both bars are labelled with their value and the vertical axis starts at zero.

Exercise 7 — Visualise — your turn

Your turn — no scaffold. A reader accepts that the nonmetro rate is higher and asks the obvious follow-up: what kind of firearm death is driving it?

Build a table called share giving, for each metro group, its total firearm deaths, its total firearm suicides, and suicides as a percentage of firearm deaths. Treat a suppressed suicide count as zero. Then draw a figure that lets a reader compare the two percentages, and print the table.

Requirements the figure must meet

A labelled vertical axis · a title · a baseline at zero · each bar labelled with its value · colours drawn from the course palette (#003594, #8a6400).

# Build the table first, the figure second.
#
# For the table: a suppressed suicide count is missing, not zero, so
# .fillna(0) converts it before you total it. Then group `joined` by metro
# and .agg() the deaths and the suicides, exactly as in Exercise 5. The
# percentage is suicides / deaths * 100, added as a new column.
# Name that column suicide_share_pct.
#
# For the figure: fig, ax = plt.subplots(...), then ax.bar() with the two
# category names and the two values. ax.bar_label(bars, fmt="%.1f") writes
# the value above each bar. ax.set_ylim(0, 80) fixes the baseline at zero —
# starting a bar axis above zero makes a small difference look enormous,
# which Module 5 called axis truncation.
#
# Keep a reference to the axes in `ax`; the check inspects it.

joined["suicides"] = joined["firearm_suicide"].fillna(0) share = joined.groupby("metro").agg(deaths=("firearm_deaths", "sum"), suicides=("suicides", "sum")) share["suicide_share_pct"] = share["suicides"] / share["deaths"] * 100 fig, ax = plt.subplots(figsize=(6, 3.2)) bars = ax.bar(["Metro", "Nonmetro"], [share.loc["Metro", "suicide_share_pct"], share.loc["Nonmetro", "suicide_share_pct"]], color=["#003594", "#8a6400"], edgecolor="white") ax.bar_label(bars, fmt="%.1f", padding=3) ax.set_ylim(0, 80) ax.set_ylabel("Suicide share of firearm deaths (%)") ax.set_title("Suicide share of 2023 firearm deaths, 42 reporting PA counties") plt.tight_layout() plt.show() print(share.round(1))

joined["suicides"] = joined["firearm_suicide"].fillna(0)
share = joined.groupby("metro").agg(deaths=("firearm_deaths", "sum"),
                                    suicides=("suicides", "sum"))
share["suicide_share_pct"] = share["suicides"] / share["deaths"] * 100

fig, ax = plt.subplots(figsize=(6, 3.2))
bars = ax.bar(["Metro", "Nonmetro"],
              [share.loc["Metro", "suicide_share_pct"],
               share.loc["Nonmetro", "suicide_share_pct"]],
              color=["#003594", "#8a6400"], edgecolor="white")
ax.bar_label(bars, fmt="%.1f", padding=3)
ax.set_ylim(0, 80)
ax.set_ylabel("Suicide share of firearm deaths (%)")
ax.set_title("Suicide share of 2023 firearm deaths, 42 reporting PA counties")
plt.tight_layout()
plt.show()

print(share.round(1))

Exercise 8 — Describe the figure

Exercise 8 — Describe the figure

Worked example. A figure that cannot be read aloud is unavailable to part of its audience and to anyone whose image fails to load. A description states the chart type, the axes, and the values — not the conclusion.

Exercise 8 — Describe the figure — your turn

Your turn — no scaffold. Someone using a screen reader will hear your description instead of seeing the figure you drew in Exercise 7. Write it into a string called share_description, and print it.

They must be able to reconstruct the figure from your words alone: what kind of chart it is, what the two groups are, what quantity is being shown, both of its values, and where the vertical axis begins. Describe what is on the page, not what it means — the conclusion is the reader’s to draw.

# Two or three sentences, written for someone who cannot see the screen.
#
# Work through what they would need: the chart type ("bar chart"), the two
# categories (metro and nonmetro), the quantity (suicide share of firearm
# deaths, as a percentage), both values (62.1 and 54.7), and the fact that
# the axis starts at zero — that last one tells them the comparison has not
# been visually exaggerated.
#
# Resist writing the conclusion. "Nonmetro counties have a higher suicide
# share" is your interpretation; the reader is entitled to reach it, or
# not, from the same numbers a sighted reader gets.

share_description = ( "Bar chart of the suicide share of 2023 firearm deaths in the 42 reporting " "Pennsylvania counties. The metro bar reaches 54.7 percent and the nonmetro " "bar reaches 62.1 percent. The vertical axis starts at zero." ) print(share_description)

share_description = (
    "Bar chart of the suicide share of 2023 firearm deaths in the 42 reporting "
    "Pennsylvania counties. The metro bar reaches 54.7 percent and the nonmetro "
    "bar reaches 62.1 percent. The vertical axis starts at zero."
)

print(share_description)

Where this goes

In an ordinary Quarto document, this text is the figure’s alt text. In these slides it is printed under the figure as visible text, because the slide format does not carry alt text through to the browser.

Exercise 9 — Report the finding

Exercise 9 — Report the finding

Worked example. A reported number is assembled from the pipeline’s own variables rather than typed by hand. A number typed by hand stops being a result the moment the data changes.

Exercise 9 — Report the finding — your turn

Your turn — no scaffold. Write the paragraph a reader would see, into a string called finding, and print it.

It must state, in prose: how many counties reported a count and out of how many; both pooled rates, identified by group; how many deaths the nonmetropolitan estimate rests on and across how many counties; and that counties reporting no total were excluded from both the numerator and the denominator.

Every number must be computed, not typed

Build the sentence with an f-string reading from summary and reporting. A number typed by hand stops being a result the moment the data changes.

# Pull each number out of the pipeline first, give it a name, then write
# one f-string that reads those names. Something like:
#
#   n_reporting   = len(reporting)
#   nonmetro_rate = summary.loc["Nonmetro", "rate_per_100k"]
#
# and so on for the metro rate, the nonmetro death count, and the nonmetro
# county count. Inside an f-string, {nonmetro_rate:.1f} prints one decimal.
#
# Report the count behind the rate, not just the rate. A reader can then
# judge how much the estimate would move if one death had landed in a
# different year. That sentence is what this course asks for instead of a
# p-value.
#
# The test of whether you did this right: change the year at the top of
# your pipeline to 2022 and run everything again. A sentence built from
# f-strings updates itself. A sentence with the numbers typed in keeps
# claiming 18.5 while the table above it says something else.

n_reporting = len(reporting) n_total = len(f2023) nonmetro_rate = summary.loc["Nonmetro", "rate_per_100k"] metro_rate = summary.loc["Metro", "rate_per_100k"] nonmetro_deaths = int(summary.loc["Nonmetro", "deaths"]) nonmetro_counties = int(summary.loc["Nonmetro", "counties"]) finding = ( f"Among the {n_reporting} of {n_total} Pennsylvania counties reporting a 2023 " f"firearm death count, the pooled rate was {nonmetro_rate:.1f} per 100,000 in " f"nonmetropolitan counties against {metro_rate:.1f} in metropolitan counties. " f"The nonmetropolitan estimate rests on {nonmetro_deaths} deaths across " f"{nonmetro_counties} reporting counties. Counties reporting no total are " f"excluded from both the numerator and the denominator." ) print(finding)

n_reporting = len(reporting)
n_total = len(f2023)
nonmetro_rate = summary.loc["Nonmetro", "rate_per_100k"]
metro_rate = summary.loc["Metro", "rate_per_100k"]
nonmetro_deaths = int(summary.loc["Nonmetro", "deaths"])
nonmetro_counties = int(summary.loc["Nonmetro", "counties"])

finding = (
    f"Among the {n_reporting} of {n_total} Pennsylvania counties reporting a 2023 "
    f"firearm death count, the pooled rate was {nonmetro_rate:.1f} per 100,000 in "
    f"nonmetropolitan counties against {metro_rate:.1f} in metropolitan counties. "
    f"The nonmetropolitan estimate rests on {nonmetro_deaths} deaths across "
    f"{nonmetro_counties} reporting counties. Counties reporting no total are "
    f"excluded from both the numerator and the denominator."
)

print(finding)

Restart and run all

Before this notebook counts as finished, restart the runtime and run every cell from the top. A result that only exists in memory is not a result.

Session outcomes

1. The seven-stage pipeline — import, clean, join, group, summarise, visualise, report

2. Suppression as a protective mechanism, and its non-random pattern

3. The silent denominator error, and the fact that it reverses this finding

4. Notebook hygiene: imports at the top, relative paths, seeds, restart-and-run-all

5. Uncertainty stated verbally — counts, coverage, and consistency across years

6. The five-move report, and computer output integrated into sentences

The finding, paid off

Nonmetro 18.5 per 100,000 against metro 14.1, with suicides at 62.1% of nonmetro firearm deaths against 54.7% in metro counties. Handled carelessly, the same data reports 8.2 against 13.6 and confirms the opposite.

The final project applies this pipeline to a question of the student’s own. The seven stages do not change — only the question and the data do.

Close by paying the promise off out loud, because it was made seventy-five minutes ago and they earned it.

Everyone in this room believed, at the start, that firearm death is an urban problem. Pennsylvania’s data says eighteen point five against fourteen point one, the other direction, and the suicide share tells you why. That is the only finding this course has shown them that contradicts what they walked in with — every other one confirmed what a public health student already suspects. Say that explicitly. It is why this one was saved for last.

And then the harder half, which is the real lesson. The wrong version of that analysis produced eight point two against thirteen point six. Clean code, no errors, a table that looks exactly as trustworthy as the right one — and it confirms the stereotype they arrived with. That is the shape of the most dangerous kind of mistake: it runs, and it agrees with you.

For the project, the sentence to leave them with is that the pipeline does not change. New question, new file, same seven stages. They have the whole skeleton now.

Practical close. Next session is the eleventh of November nineteenth, Final Project Advancement — bring a question and the file they intend to use, not a finished analysis. Presentations are December first and third. And the writeup is thirty percent, which is worth saying one final time.