Module 3: NumPy

PUBHLT 0411 · Python for Public Health Data Analysis

Soumik Purkayastha

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

October 22, 2026

Nine micrograms.

Sixty-seven counties and one federal threshold. No loops allowed!

The EPA annual standard for PM2.5 is 9 micrograms per cubic metre

Twenty-two of Pennsylvania’s 67 counties exceed it. One exceeds it by a wide margin.

The question this session opens

Ranking 67 counties against a threshold, and measuring how far each sits from the state’s typical value, requires arithmetic on every county at once.

This session establishes the array: one object holding all 67 values, arithmetic that applies to every element in a single expression, and selection by condition rather than by position.

Required from Module 2: for, the accumulator pattern, def and return, and the distinction between counts and rates.

Most counties sit near the standard; Allegheny sits far above it

Figure 1. PM2.5 annual averages for all 67 Pennsylvania counties, drawn as a violin with a box plot inside and one dot per county below. The median is 8.5 and the middle half runs from 7.5 to 9.2, so the bulk of the state straddles the dashed EPA line at 9.0. Twenty-two counties, in blue, sit above it. Allegheny, at 14.1, sits alone at the right edge.

The violin shows the shape; the box shows the five-number summary; the dots show that every shape is built from 67 real counties. Point at each layer in that order.

The story is two-part. First, the standard cuts through the middle of the state, not the tail — a third of the counties are over. Second, Allegheny is a different kind of over: the next highest county is Lancaster at 11.1, three full units below.

Producing this figure needs the median, the quartiles, and a count of values above 9. By the end of today each of those is one line of NumPy.

A NumPy array holds many values of one type

1. np.array() builds an array from a list

2. .dtype reports the single type every element shares

3. .shape reports the size along each dimension, as a tuple

import numpy as np is the universal convention

Every NumPy example in every textbook and every answer online uses the name np. The alias is not required by the language, and departing from it makes code harder for others to read.

These five values are real: Adams, Allegheny, Armstrong, Beaver, and Bedford, in the order the file lists them. Allegheny is the 14.1.

.shape returning (5,) with the trailing comma catches people. That is Python’s notation for a one-element tuple, and it is distinguishing (5,) the tuple from (5) the number five in parentheses. Say it once; it recurs constantly in Module 4.

The one-type rule is the trade. A list holds anything, an array holds one type, and that restriction is exactly what buys the speed on the next slide.

An array applies arithmetic to every element at once

A list requires a loop. An array does not.

Each operation produces a new array of the same shape, with the operation applied element by element. This is vectorisation.

The same operation on a list does something else entirely

[9.1, 14.1] * 2 repeats the list, producing four elements. np.array([9.1, 14.1]) * 2 doubles each value, producing two. The operator is the same; the type decides the meaning.

Reading a real file into an array

np.genfromtxt() reads a CSV directly into arrays. Column 0 holds names; column 1 holds the measurement.

skip_header=1 discards the header row. usecols= selects a column. Text needs dtype=str; numbers are read as float64 by default.

Source

County Health Rankings, PM2.5 annual average by county, all 67 Pennsylvania counties.

First real file of the course. Worth marking as a milestone — from here on the numbers on screen are read from disk, not typed in.

Two calls for two columns is clumsy and they should notice. That clumsiness is the argument for Module 4: pandas reads the whole table in one call and keeps the columns together. Say that explicitly, because “this is annoying” is the correct reaction and it is about to be solved.

The parallel-arrays arrangement — counties[3] and pm25[3] describing the same county — is the thing to watch. Nothing enforces the alignment. Sort one without the other and the data is silently wrong, which is a real bug and a preview of why a DataFrame exists.

Indexing and slicing follow Module 1’s rules

Zero-indexing, exclusive endpoints, negative indices, and strides are identical to string and list slicing.

Nothing new here, and that is the point worth stating: the slicing rules learned on strings in Module 1 transfer without modification. One rule, three data structures.

pm25[0] is Adams at 9.1, pm25[1] is Allegheny at 14.1, and pm25[-1] is York at 8.9, because the file is alphabetical.

The stride [::10] takes every tenth county. It has no analytical meaning here — it is a syntax demonstration, so do not invent a reason for it.

A two-dimensional array is indexed [row, column]

rates[1, 3] selects one element. rates[2] selects a whole row. rates[:, 4] selects a whole column — the colon means “every row.”

Rows and columns carry meaning

Here a row is a county and a column is a year. The shape (4, 5) states that directly: four counties, five years.

Real numbers from the overdose file, and the same four counties that have appeared since Module 1. Philadelphia’s 91.4 in 2022 is the highest cell in the table.

rates[:, 4] is the line to slow down on. The colon in the row position means all rows, and the 4 in the column position means the fifth year — so it reads out 2023 for every county. Read it aloud as “every row, column four.”

Centre County’s row is the one to point at. It runs 5.1, 8.9, 7.6, 8.9, 4.4 while Philadelphia runs in the seventies and eighties. That is the Module 1 comparison, now sitting inside one object.

axis= chooses the direction a summary collapses

1. No axis= collapses everything to a single number

2. axis=1 collapses the columns, leaving one value per row — a per-county mean

3. axis=0 collapses the rows, leaving one value per column — a per-year mean

axis= names the dimension removed, not the one kept

axis=0 runs down the rows and returns a value per column. Reading it as “axis 0 means rows, so I get row answers” is the standard mistake, and it produces a plausible array of the wrong length.

This is the hardest idea in the module and it deserves the time.

The reliable way to check is the shape. Start with (4, 5). axis=1 removes the 5, leaving 4 numbers — one per county. axis=0 removes the 4, leaving 5 numbers — one per year. Tell them to check the length of what comes back whenever they are unsure, because the wrong axis usually still returns an array and stays silent.

The per-county means are 54.4, 82.3, 42.5, and 7.0. The per-year means climb from 41.9 in 2019 to 50.0 in 2022 and then fall slightly. Four counties is too few to call that a trend, and saying so is the honest move.

This is the same axis= that .groupby() generalizes in Module 4.

Summary statistics are methods on the array

.mean(), .std(), .min(), and .max() are methods on the array. np.median() and np.percentile() are functions taking the array as an argument.

.std() divides by n, not by n − 1

NumPy’s default is the population standard deviation. R’s sd() and most statistics courses use the sample version, which is pm25.std(ddof=1). On 67 counties the two differ by under one percent, but the default is not the one a statistics course teaches.

The method-versus-function split is arbitrary and students will ask why. The honest answer is history, not design — some operations became array methods and some stayed as module functions. Do not manufacture a principle. When in doubt, try the method and fall back to np.something().

The ddof point is worth the thirty seconds because this course runs alongside a statistics sequence. Default ddof=0 gives 1.28 here, ddof=1 gives 1.29. The numbers on the payoff slide use the default, which is what .std() produces.

Mean 8.49 against median 8.50 — nearly identical, so this distribution is close to symmetric apart from Allegheny. That is worth contrasting with the overdose rates in the lab, where the mean sits well above the median.

A boolean mask selects elements by condition

A comparison on an array produces an array of booleans. Indexing with that array returns the elements where it is True.

.sum() on a boolean array counts the True values, because True equals 1.

The mask replaces the counting loop

Module 2 counted with a loop, an if, and a counter initialized outside. The mask states the condition once and the count follows from it.

Point back at Module 2’s counting exercise explicitly. Five lines of loop become one comparison and a .sum(). That is the payoff for the whole module and it is best shown as a before-and-after.

The True-is-1 trick is worth a sentence. It is not a NumPy convenience, it is how Python booleans work — True + True is 2 at the interpreter prompt. Try it.

.mean() on a boolean array gives the proportion rather than the count, which is often what is actually wanted. Mention it; the lab uses it.

Masking the counties above the standard

The same mask indexes both arrays, because the two are aligned element by element. ~ negates a mask.

Masks use &, |, and ~ — not and, or, not

and and or require a single true-or-false value and raise ValueError on an array. Each comparison also needs its own parentheses: (pm25 > 9.0) & (pm25 < 11.0).

Twenty-two counties, 32.8% of the state. The eight printed are alphabetical, not ranked — Adams, Allegheny, Beaver, Berks, Blair, Bucks, Butler, Cambria.

counties[over] is the move that makes masks useful. The mask was computed from the measurements and applied to the names, and it works only because position 12 means the same county in both arrays. That alignment is a promise nothing in the code enforces.

The &-versus-and warning is on the slide because it is the single most common NumPy error and it comes back in Module 4 with pandas, where the parentheses trip people even harder. It is worth demonstrating the ValueError so the message is familiar: “The truth value of an array with more than one element is ambiguous.”

The two conditional means, 9.80 against 7.85, are the sentence a report would carry.

np.where() builds a new array from a condition

np.where(condition, value_if_true, value_if_false) returns a full-length array, choosing element by element.

A mask filters; np.where() labels

pm25[pm25 > 9.0] returns 22 values. np.where(pm25 > 9.0, ...) returns 67, one per county. Selecting and labelling are different jobs.

The filter-versus-label distinction is the one to land. Both start from the same comparison and they answer different questions — “which counties” versus “what does each county get called.”

Nesting works for three or more bands: np.where(a, "High", np.where(b, "Middle", "Low")). It is legal, it is common, and past two levels it becomes unreadable. Show it once, then say pandas has a better tool for it in Module 4.

np.where() with one argument does something completely different — it returns the indices where the condition holds. Do not open that door unless someone asks.

Standardizing expresses each value in standard deviations

A z-score measures distance from the mean in units of the standard deviation.

The whole transformation is one line, because subtraction and division apply to every element.

A z-score is a statement about this collection

The value describes a county’s position among Pennsylvania counties. It carries no information about health effects, and a z-score near zero can still exceed a regulatory standard.

Two counties clear +2: Allegheny at 4.37 and Lancaster at 2.03. Cumberland and Delaware are next, both at 1.80, and both miss the cut. The mean of z prints as effectively zero and the standard deviation prints as 1.00 — that is the arithmetic working, and it is a good self-check for them to remember.

The caveat on the slide matters here. Philadelphia’s PM2.5 is 8.7, which is a z of +0.16 — almost exactly typical for Pennsylvania — and it is still under the federal standard. Two different questions: how does this county compare to its neighbours, and does it meet the rule. Standardizing answers the first only.

Do not extend this into an exposure-disparity claim. County-level averages cannot support one, and in Pennsylvania the correlation with poverty runs the wrong way for the story students expect.

np.nan marks a value that is absent

Any arithmetic touching nan produces nan. One missing value makes the whole mean missing, which is the safe behaviour: the result announces the gap rather than hiding it.

nan is never equal to anything, including itself

np.nan == np.nan evaluates to False. Testing for missingness requires np.isnan().

The suppression is real and worth explaining properly. Pennsylvania withholds county counts below ten because a rate computed from three deaths is unstable and because small counts can identify people. Seven of the 67 counties are suppressed for 2023.

nan propagating is a feature, not a bug. A mean that silently ignored the gap would report a state figure that quietly excluded seven rural counties. Getting nan back is the array saying “you have not decided what to do about this yet.”

The self-inequality is genuinely strange and they should just memorize it. It follows from the floating-point standard: nan means “not a number,” and no comparison with an unknown quantity can be true.

This is the seed of Module 4’s whole missing-data section, and it is the same judgment as Module 2’s .get(key, 0). Missing is not zero.

The nan-aware functions skip the gaps

np.nanmean(), np.nanmedian(), np.nansum(), and np.nanstd() compute over the present values.

Skipping a value is a decision, and it belongs in the writeup

np.nanmean() computes the mean of the counties that reported. Where suppression is concentrated in small rural counties, that mean describes a different population than the one requested, and the count of omitted counties is part of the result.

Report the denominator. Always. “30.1 per 100,000 across the 60 counties that reported” is honest; “30.1 per 100,000” alone is not, and the difference costs nothing to write.

Across all 67 counties for 2023 the real figures are: 60 reporting, mean 30.1, median 26.8. The mean above the median tells you the distribution has a long right tail — a few counties far worse than typical — which is exactly the opposite of the PM2.5 shape from earlier.

The instinct students have is to fill the gaps with zero. Ask what that would assert: that nobody in Cameron County died of an overdose in 2023. The state suppressed the number precisely because it is not zero. Module 4 has the full treatment.

A slice is a view; modifying it modifies the original

A slice of an array does not copy. It is a view onto the same memory, and writing through it changes the original array.

Slicing a list copies; slicing an array does not

list[:3] produces an independent copy. array[:3] produces a view. The syntax is identical and the behaviour is not.

Resolution: .copy() makes an independent array — first_three = pm25[:3].copy().

This is the module’s real bug source and the reason it is on the spec. The failure is silent: a value changes in an array the code never appears to touch, and the wrong number propagates through everything downstream.

Run it with .copy() added and show that pm25 keeps its 9.1. One character of difference, completely different behaviour.

Why NumPy does this: arrays are meant to be large, and copying a slice of a million-element array on every operation would be ruinous. The view is a performance decision, and the copy is available whenever it is wanted.

One useful asymmetry worth knowing: boolean-mask indexing does copy. pm25[pm25 > 9] returns new memory, so writing to it is safe. Slices are views, masks are copies.

The outlier.

Putting the session together.

Exercise: standardize and rank

# The blank is mean. A z-score subtracts the mean, then divides by the
# standard deviation. Subtracting the median would centre the values but
# would not produce a z-score.
#
# .argmax() returns the *position* of the largest value, which then indexes
# both arrays — the same alignment used throughout this session.

import numpy as np counties = np.genfromtxt( "data/pa_air_pollution.csv", delimiter=",", skip_header=1, usecols=0, dtype=str ) pm25 = np.genfromtxt( "data/pa_air_pollution.csv", delimiter=",", skip_header=1, usecols=1 ) z = (pm25 - pm25.mean()) / pm25.std() worst = z.argmax() print(f"Highest PM2.5: {counties[worst]}") print(f" concentration: {pm25[worst]} ug/m3") print(f" z-score: {z[worst]:+.2f}") print(f"Counties over the 9.0 standard: {(pm25 > 9.0).sum()}")

import numpy as np

counties = np.genfromtxt(
    "data/pa_air_pollution.csv", delimiter=",", skip_header=1, usecols=0, dtype=str
)
pm25 = np.genfromtxt(
    "data/pa_air_pollution.csv", delimiter=",", skip_header=1, usecols=1
)

z = (pm25 - pm25.mean()) / pm25.std()

worst = z.argmax()

print(f"Highest PM2.5: {counties[worst]}")
print(f"  concentration: {pm25[worst]} ug/m3")
print(f"  z-score:       {z[worst]:+.2f}")
print(f"Counties over the 9.0 standard: {(pm25 > 9.0).sum()}")

Fill this in with the room. It is the whole module in twelve lines — read a real file, vectorized arithmetic, a summary method, and a mask.

.argmax() is new on this slide and it is worth naming. .max() gives the largest value, .argmax() gives its position, and the position is what indexes the names array. That distinction comes back constantly.

The :+.2f format forces the plus sign. Small presentation detail, and presentation is graded.

Ask for the county before running. They will say Allegheny and they will be right. The magnitude is the part nobody predicts.

Allegheny County

Figure 2. Each dot is one Pennsylvania county’s annual average PM2.5. Sixty-six counties fall between 5.5 and 11.1; Allegheny sits alone at 14.1, a blue diamond well to the right of every other county. Twenty-two counties exceed the EPA standard of 9.0, marked in gold, and all but Allegheny cluster between 9.1 and 11.1. Two arrows compare the gaps: 3.0 from Lancaster up to Allegheny, and 2.6 from the state mean of 8.49 up to Lancaster.

Allegheny County’s PM2.5 average is 14.1 micrograms per cubic metre — a z-score of +4.37.

The state mean is 8.49. The next-highest county, Lancaster, is 11.1.

What the magnitude means

A z-score above 4 places Allegheny more than four standard deviations from the typical Pennsylvania county. The gap between Allegheny and second place is larger than the gap between second place and the state mean.

Twenty-two counties exceed the 9.0 standard. Twenty-one of them cluster between 9.1 and 11.1. Allegheny sits alone.

Land it and stop.

The distribution matters more than the single number. Strip Allegheny out and Pennsylvania’s counties run from 5.5 to 11.1 — a normal-looking spread. Allegheny is not the top of a gradient, it is a separate case, and that is what a z-score above 4 is telling you.

Say plainly what this is and is not. It is a county-wide annual average against a federal standard. It is not a statement about any neighbourhood, and PM2.5 varies enormously within Allegheny County — the Mon Valley is not Sewickley. County averages are the resolution this data has, and the honest move is to name that limit rather than to over-read the number.

Then close the loop on Module 2: this took one line of arithmetic across 67 counties. The equivalent loop would be fifteen lines, and it would be slower.

Challenge: the outlier hides the others

Allegheny inflates the standard deviation for every county. Remove it and recompute.

# The first blank builds a mask that is True for every county EXCEPT
# Allegheny: the comparison operator for "not equal". The second blank is
# the same summary method used in the z-score above.
#
# counties[keep] applies the same mask to the names, so position `top`
# in z_rest lines up with position `top` in counties[keep].

keep = counties != "Allegheny" rest = pm25[keep] z_rest = (rest - rest.mean()) / rest.std() top = z_rest.argmax() print(f"Without Allegheny, highest: {counties[keep][top]}, z = {z_rest[top]:+.2f}") print(f"Beyond 2 SD, all 67: {((z > 2) | (z < -2)).sum()}") print(f"Beyond 2 SD, without it: {((z_rest > 2) | (z_rest < -2)).sum()}")

keep = counties != "Allegheny"
rest = pm25[keep]
z_rest = (rest - rest.mean()) / rest.std()

top = z_rest.argmax()
print(f"Without Allegheny, highest: {counties[keep][top]}, z = {z_rest[top]:+.2f}")
print(f"Beyond 2 SD, all 67:        {((z > 2) | (z < -2)).sum()}")
print(f"Beyond 2 SD, without it:    {((z_rest > 2) | (z_rest < -2)).sum()}")

Your turn!

Lancaster’s PM2.5 did not change, yet its z-score rose from +2.03 to +2.47. Why? And why does removing one extreme county produce more counties beyond 2 SD, not fewer?

The answer: the standard deviation drops from 1.28 to 1.09 once Allegheny is gone, so every other county’s distance is now measured in a smaller unit. Lancaster moves from +2.03 to +2.47, and Cumberland and Delaware cross the line at +2.20. Three counties beyond 2 SD becomes four.

The lesson to land: a z-score is relative to the collection, as the earlier callout said. One extreme value stretches the ruler and makes everything else look ordinary. This is why an outlier is worth checking before it is used to standardise anything.

!= on an array of strings returns a mask, the same as > on numbers. They have seen ~; counties != "Allegheny" is the direct form.

Session outcomes

1. Building arrays with np.array(), and reading real ones with np.genfromtxt()

2. Inspecting .dtype, .shape, and .size

3. Indexing and slicing in one dimension, and [row, column] in two

4. Vectorised arithmetic — one expression applied to every element

5. Summarising with .mean(), .std(), np.median(), and np.percentile()

6. Collapsing a chosen dimension with axis=

7. Selecting with a boolean mask, and labelling with np.where()

8. Standardising to z-scores

9. Handling np.nan with np.isnan() and the nan-aware functions

10. Distinguishing a view from a copy

Preparation for the lab, Tuesday 10/27

The lab reads both course files, masks against the standard, standardises overdose rates, and handles the suppressed counties.

Module 4 attaches these arrays to labelled columns and calls the result a DataFrame.

The single sentence for the whole module: an array applies one operation to every element, so the loop disappears.

Close on Module 4 and be concrete about what it fixes. Today the county names and the measurements sat in two arrays held in alignment by nothing but care. A DataFrame keeps them in one object, so sorting or filtering moves them together and the alignment cannot break. Every complaint about the two genfromtxt calls is answered on Thursday.

That’s all folks!

See you on Tuesday!