Data Science

A cricket commentator who says “this batter scores most of his runs on the off side” is doing a little data science. So is a farmer who keeps a notebook of when she planted and how much she harvested, and a mobile company that looks at when people use the most data to decide where to build towers. Each starts with a question, looks at data, and draws a conclusion.

In Grade 7 you met the steps of a data investigation (Figure 4.1). They are the same in data science, and the same whether the data is 30 answers from your class or millions of rows from a whole country:

What changes with a computer is the size of what you can do, and one new tool. A spreadsheet is good for a few hundred rows that you can look at. A program can clean a million rows the same way every time, and you can read the program to see exactly what was done, which you can’t do with a spreadsheet whose cells were changed by hand. And when the question is about something that hasn’t happened yet, or can’t be tried for real (“what if the hospital opened a second counter?”), a program can make up data that behaves like the real thing and try it out: a simulation.

A person who does this work as a job is called a data scientist or a data analyst. They work in banks, mobile companies, hospitals, government offices, newspapers and research centers, in Nepal and everywhere else.

In Grade 9 you read a CSV file into a list of dictionaries, one dictionary per row, and wrote loops to add up and compare. That works, but every question needs a new loop. pandas is a library that does it for you. It keeps a whole table in one variable, called a DataFrame: a table with named columns and numbered rows, like a sheet in a spreadsheet. Each column is a Series, a list of values of one type with its own name.

Unit 1 promised that some libraries have a short name everyone uses. pandas is one: every program, every book and every answer online imports it as pd:

import pandas as pd

Reading a table

Our first data is the table you used in Grade 9’s machine learning unit: Nepal’s 77 districts, one row each, with each district’s belt (Mountain, Hill or Terai), its height (the middle height of its land, in meters), its slope (how steep its land is on average, in degrees), its rain (the rain in an average year, in millimeters) and its people (people per square kilometer, from the 2021 census).

pd.read_csv reads the whole file into a DataFrame in one line. head() shows its first five rows, and shape says how big it is:

import pandas as pd districts = pd.read_csv("districts.csv")print(districts.head())print(districts.shape)

       district      belt  height  slope  rain  people0        Achham      Hill    1308     26  1987     1351  Arghakhanchi      Hill     990     21  2324     1432       Baglung      Hill    2235     28  2730     1363       Baitadi      Hill    1600     26  2245     1624       Bajhang  Mountain    3261     31  2272      55(77, 6)

pandas has numbered the rows from 0 on the left: that is the table’s index. (77, 6) means 77 rows and 6 columns. The first time you press Run, the page has to download pandas, so it takes a little longer.

Picking columns and rows

A column is picked by its name in square brackets, as with a dictionary. A column of numbers has methods that work on every value at once, with no loop:

import pandas as pd districts = pd.read_csv("districts.csv")rain = districts["rain"]print("Average rain:", round(rain.mean()), "mm")print("Least:", rain.min(), "mm, most:", rain.max(), "mm") wettest = districts.sort_values("rain", ascending=False)print(wettest.head(3)) terai = districts[districts["belt"] == "Terai"]print(len(terai), "districts are in the Terai")

Average rain: 2198 mmLeast: 792 mm, most: 3777 mm         district      belt  height  slope  rain  people63  Sankhuwasabha  Mountain    2625     29  3777      4633          Kaski      Hill    2131     27  3630     28838        Lamjung      Hill    1924     28  3479      9421 districts are in the Terai

Two lines need a closer look.

• sort_values("rain", ascending=False) gives a new DataFrame with the rows in order of rain, the most first. The original, districts, is left as it was.
• districts["belt"] == "Terai" compares every value in the column with "Terai" and gives a Series of True and False, one per row. Putting that Series inside districts[...] keeps only the rows where it is True. This is how pandas filters rows, like the filter in a spreadsheet or SQL’s WHERE in Unit 3. Conditions join with & (and) and | (or), each in its own brackets: districts[(districts["belt"] == "Hill") & (districts["rain"] > 3000)].

Groups, the middle and the spread

“Which belt is most crowded?” is a question about groups: split the rows by belt, then work out one number for each group. groupby does exactly that (it is SQL’s GROUP BY, from Unit 3):

import pandas as pd districts = pd.read_csv("districts.csv")by_belt = districts.groupby("belt")["people"]print(by_belt.agg(["count", "mean", "median"]).round())

          count   mean  medianbelt                          Hill         40  384.0   146.0Mountain     16   42.0    40.0Terai        21  522.0   551.0

Look at the Hill row. The mean (the average) is 384 people per km², but the median, the middle district once they are in order, is only 146. Which is the better summary of a hill district? Most hill districts (29 of the 40) have between 100 and 250 people per km². Three districts in the Kathmandu valley have over 1,000, Kathmandu itself nearly 5,000, and those three pull the mean far up. A value far from the rest is called an outlier. The median doesn’t care how far out an outlier is, only which side of the middle it falls on, so with outliers the median says more about a typical member of the group. Unit 1 promised you this: here is where the median earns its place.

The spread matters too. describe() gives the main summaries of every number column at once:

import pandas as pd districts = pd.read_csv("districts.csv")print(districts.describe().round())

       height  slope    rain  peoplecount    77.0   77.0    77.0    77.0mean   1570.0   20.0  2198.0   351.0std    1342.0   10.0   543.0   695.0min      93.0    1.0   792.0     2.025%     340.0   10.0  1869.0    93.050%    1308.0   23.0  2060.0   166.075%    2235.0   28.0  2561.0   318.0max    4933.0   32.0  3777.0  4936.0

std is the standard deviation you met in Unit 1: how far, roughly, the values are from their mean. 25%, 50% and 75% are the quartiles: a quarter of the districts are below the 25% value, half below the 50% value (the median), three quarters below the 75% value. So half of all districts get between 1,869 and 2,561 mm of rain a year. Two districts with the same mean can still be very different places: the mean tells you the middle, the spread how much to expect around it.

Real data almost never arrives clean. You saw the usual problems in Grade 7: the same answer in different spellings, numbers written as words, a wrong unit, missing values, repeats and impossible values. In a spreadsheet you fixed them by hand. With a program you write each fix once, as a step, and the steps can be read, checked and run again on next year’s data. Cleaning data with a program is often most of the work of a data scientist.

Here is a task where it matters. You want to know how many people of Bagmati Province live in each belt. districts.csv knows each district’s belt, but not its population. A report gives the population of Bagmati’s 13 districts from the 2021 census, and you copy its table into a file. The numbers are the census’s; the mess is the kind that copying from reports really makes (we put it in for practice):

district,population Kathmandu,"20,41,587"lalitpur,5.52 lakhBhaktapur,"4,32,132"Kavre,"3,64,039"Sindhupalchowk,"2,62,624"Dolakha,"1,72,767"Ramechhap,"1,70,302"Sindhuli,"3,00,026"Makawanpur,"4,66,073"Chitwan,"7,19,859"Dhading ,"3,25,710"Nuwakot,"2,63,391"Rasuwa,-Bhaktapur,"4,32,132"

(The commas are in the Nepali places, 20,41,587, twenty lakh forty-one thousand five hundred and eighty-seven. A CSV file puts a value with commas in quotes, so that its commas aren’t taken for the commas between columns.)

Look before you clean

The first step of cleaning is to look, and to let pandas tell you what it sees. dtypes gives each column’s type:

import pandas as pd people = pd.read_csv("bagmati_population.csv")print(people.dtypes)print(people["district"].tolist())print(people["population"].sum())

district      strpopulation    strdtype: object[' Kathmandu', 'lalitpur', 'Bhaktapur', 'Kavre', 'Sindhupalchowk', 'Dolakha', 'Ramechhap', 'Sindhuli', 'Makawanpur', 'Chitwan', 'Dhading ', 'Nuwakot', 'Rasuwa', 'Bhaktapur']20,41,5875.52 lakh4,32,1323,64,0392,62,6241,72,7671,70,3023,00,0264,66,0737,19,8593,25,7102,63,391-4,32,132

Three warnings in four lines:

• The population column’s type is str, text, not a number. pandas makes a column a number only when every value in it reads as a number, and "20,41,587", "5.52 lakh" and "-" don’t.
• The list of names shows what the table hides: a space before Kathmandu and after Dhading, a small l in lalitpur, and Bhaktapur twice.
• sum() of text doesn’t fail. It joins the strings end to end, as + does for strings, into one long meaningless string. If you hadn’t looked, a program printing a “total” would have printed nonsense.

Cleaning step by step

Here are the four cleaning steps, written as one function in a module, the way Unit 1 built programs from modules. Each step is one or two lines, with a comment saying why:

"""Cleaning the population table copied from the report.""" import pandas as pd # other spellings of district names, and the spelling districts.csv usesSPELLINGS = {"Kavre": "Kavrepalanchok",             "Sindhupalchowk": "Sindhupalchok",             "Makawanpur": "Makwanpur"}  def to_number(text):    """A population written as text, as a whole number, or None if missing."""    text = text.strip()    if text == "-":        return None                      # missing: not zero    if text.endswith("lakh"):            # "5.52 lakh" is 5.52 × 100,000        return round(float(text.removesuffix("lakh")) * 100_000)    return int(text.replace(",", ""))    # "20,41,587" → 2041587  def clean_population(file_name):    people = pd.read_csv(file_name)    # 1. no spaces at the ends of names, and a capital first letter    people["district"] = people["district"].str.strip().str.title()    # 2. one spelling for each district    people["district"] = people["district"].replace(SPELLINGS)    # 3. numbers from text; Int64 is a whole-number type that can hold a gap    people["population"] = people["population"].apply(to_number).astype("Int64")    # 4. a row that is exactly the same as an earlier one is a repeat    people = people.drop_duplicates()    return people

Some new tools in it:

• .str gives a text column all of Python’s string methods (Grade 9), done to every value: .str.strip() takes spaces off both ends, .str.title() makes the first letter of each word a capital.
• replace(SPELLINGS) swaps every value that is a key in the dictionary for its value, and leaves the rest alone.
• apply(to_number) calls our own function once for each value and makes a new column of what it returns. astype("Int64") says the column holds whole numbers; pandas shows a missing one as <NA>.
• drop_duplicates() keeps only the first of rows that are the same in every column.

Now use it:

import cleaning people = cleaning.clean_population("bagmati_population.csv")print(people)print(people.dtypes)

          district  population0        Kathmandu     20415871         Lalitpur      5520002        Bhaktapur      4321323   Kavrepalanchok      3640394    Sindhupalchok      2626245          Dolakha      1727676        Ramechhap      1703027         Sindhuli      3000268        Makwanpur      4660739          Chitwan      71985910         Dhading      32571011         Nuwakot      26339112          Rasuwa        <NA>district        strpopulation    Int64dtype: object

Thirteen districts, each once, with numbers that are numbers. Rasuwa’s gap is still a gap: the report didn’t give it, and inventing a number, or writing 0, would make the totals wrong in a way nobody could see. Lalitpur’s 552,000 is only as exact as “5.52 lakh”, so it may be a few hundred out. A good cleaning program changes only what it must, and someone should write both of these down for the report.

Joining two tables

Now each district needs its belt from districts.csv. Joining two tables by a column they share is a merge in pandas: it is SQL’s JOIN … ON from Unit 3.

import pandas as pdimport cleaning people = cleaning.clean_population("bagmati_population.csv")districts = pd.read_csv("districts.csv")both = people.merge(districts, on="district", how="left")print(both[["district", "belt", "population"]])print(both.groupby("belt")["population"].sum())

          district      belt  population0        Kathmandu      Hill     20415871         Lalitpur      Hill      5520002        Bhaktapur      Hill      4321323   Kavrepalanchok      Hill      3640394    Sindhupalchok  Mountain      2626245          Dolakha  Mountain      1727676        Ramechhap      Hill      1703027         Sindhuli      Hill      3000268        Makwanpur      Hill      4660739          Chitwan     Terai      71985910         Dhading      Hill      32571011         Nuwakot      Hill      26339112          Rasuwa  Mountain        <NA>beltHill        4915260Mountain     435391Terai        719859Name: population, dtype: Int64

how="left" keeps every row of the left table, people, even one with no match, with an empty belt: Unit 3’s LEFT JOIN. Then you can check that nothing went unmatched: both["belt"].isna().sum() counts the empty belts, and it should be 0. A plain merge (how="inner", the default) would quietly drop the rows that don’t match, as a plain JOIN does.

That is why cleaning comes first. Try the merge on the uncleaned table, and " Kathmandu" with its space doesn’t equal "Kathmandu", "Kavre" doesn’t equal "Kavrepalanchok", and about four million of Bagmati’s people vanish from the totals without a single error message. Try it below.

Cleaning answers people typed

Text that people type themselves is the messiest data of all. A class answered “How do you come to school?” in a box, and got answers like these. To count them, each answer must be put into one of a few groups, by the words in it:

import pandas as pd answers = pd.Series(["Bus", "by bus ", "walk", "Walking", "on foot", "school bus",                     "Micro", "scooter with dad", "Cycle", "bicycle", "WALK", "bike", ""]) def way(text):    text = text.strip().lower()    if text == "":        return "no answer"    for word in ["walk", "foot"]:        if word in text:            return "walk"    for word in ["bus", "micro"]:        if word in text:            return "bus or micro"    for word in ["cycle", "bicycle"]:        if word in text:            return "bicycle"    return "other" print(answers.apply(way).value_counts())

bus or micro    4walk            4other           2bicycle         2no answer       1Name: count, dtype: int64

value_counts() counts how many times each value appears, most first. Every rule here is a choice that someone made: “bike” fell into “other”, because in Nepal a bike is often a motorbike; and a student who wrote “walk then bus” would count as a walker, because the walk rule comes first. So look at what fell into “other”, and write down your rules beside your results.

A table of 77 districts with six columns is too much to read. A chart shows it all at once, and lets your eyes find what a table hides. In Grade 7 you drew a scatter chart: one dot per row, placed across by one number and up by another, to see if the two go together. Each column that a chart shows is called a variable. A scatter chart shows two.

But a dot has more than a place. It has a color, a size and a shape, and each of these can show one more variable. A chart that shows three or more variables is a multivariate chart (“many variables”). With position across, position up and color, one chart shows three columns of the districts table at once.

Notice what the color added. With only height and rain, you see the hump. With the belts colored, you also see that “the Terai” and “the hills” are not one kind of place each, and that a chart of averages by belt would have hidden the hump completely.

Drawing it in Python

Matplotlib draws charts in Python. Its pyplot module is imported as plt, another name that everyone uses. Here is the chart of height and rain, colored by belt. On this page, plt.show() shows the chart under the program:

import pandas as pdimport matplotlib.pyplot as plt districts = pd.read_csv("districts.csv")colors = {"Mountain": "purple", "Hill": "green", "Terai": "orange"} for belt, group in districts.groupby("belt"):    plt.scatter(group["height"], group["rain"], color=colors[belt], label=belt) plt.xlabel("Height (m)")plt.ylabel("Rain (mm a year)")plt.title("Nepal's districts: height, rain and belt")plt.legend()plt.show()

groupby("belt") gives each belt’s name with its group of rows, so the loop draws the three belts one after another, each in its own color with its own label for the legend, the key that says which color means what. Every chart needs a title, both axes named with their units, and a legend when color means something: Grade 6’s checklist still holds.

Small multiples

A fourth variable could go in the dots’ size, but sizes are hard to compare by eye. Another way is to draw the same chart several times, once for each group, side by side, on the same scales. These are called small multiples. Here are the three belts, with slope across and people up, on a log scale:

import pandas as pdimport matplotlib.pyplot as plt districts = pd.read_csv("districts.csv")fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(9, 3.5)) for ax, belt in zip(axes, ["Mountain", "Hill", "Terai"]):    group = districts[districts["belt"] == belt]    ax.scatter(group["slope"], group["people"])    ax.set_title(belt)    ax.set_xlabel("Slope (degrees)") axes[0].set_ylabel("People per km²")axes[0].set_yscale("log")plt.show()

plt.subplots(1, 3) makes one row of three charts, axes. sharex and sharey give them the same scales, which is the whole point: a dot at the same place means the same thing in all three. Run it and compare: within the Hill belt, are flatter districts more crowded? Which three dots break the pattern?

Which chart for which question

Some questions can’t be answered from data you have. A hospital in a district town sells tickets for its outpatient department from one counter, and the line on a busy morning is long. Would a second counter help? It could open one and see, but that costs a salary for months, and some mornings are busier than others, so a few days wouldn’t tell much. Instead, you can build a simulation: a program that acts out the counter’s morning with made-up people who arrive and are served the way real people are. Then you can try one counter, two, three, on a hundred mornings each, in a second, and nobody waits in a real line.

A simulation runs on simulated data: data made up by a program, following rules that make it behave like the real thing. The rules are given as a few numbers, called parameters: here, how many people come in an hour, and how long serving each one takes.

Data made to order

Python’s random module (Unit 1) can make up numbers that follow parameters you choose. Suppose you want to test Ms. Karki’s results program from Unit 1 on a whole class, without using real students’ marks, which are private. A made-up class of 40, with marks out of 75 that average about 48 and spread about 12 either side:

import randomimport statistics random.seed(13)marks = []for i in range(40):    mark = round(random.gauss(48, 12))   # around 48, spread about 12    mark = max(0, min(75, mark))         # never below 0 or above 75    marks.append(mark) print(marks[:12])print("Average:", round(statistics.mean(marks), 1))print("Spread:", round(statistics.stdev(marks), 1))print("Lowest and highest:", min(marks), max(marks))

[47, 66, 39, 27, 51, 56, 53, 55, 47, 42, 40, 46]Average: 48.8Spread: 12.0Lowest and highest: 27 71

• random.gauss(48, 12) gives a number from the normal distribution, the bell shape that many measurements follow (heights, marks, the weights of mangoes): most near the middle, 48, fewer and fewer far from it, with a standard deviation of 12. About two thirds of its values fall within one standard deviation of the middle, here between 36 and 60.
• max(0, min(75, mark)) keeps every mark possible: a bell shape goes on for ever both ways, but marks don’t.
• random.seed(13) starts Python’s random numbers from a fixed point, so the program makes the same “random” class every time it runs, which you need to check a result or to find a bug. Change the 13 to get a different class, or leave the line out for a new one every run.

Then check that the data fits: the average and spread of the made-up marks are close to 48 and 12, but not exactly, as with any 40 real students. Simulated data that doesn’t match its parameters, or real data, is the first thing to suspect when a simulation gives a strange answer.

People arriving at a counter follow a different rule. On average 25 come in an hour, one every 2.4 minutes, but not like a clock: each person decides for themselves when to come, so sometimes three arrive within a minute and sometimes nobody comes for ten. The gaps between arrivals like that follow the exponential distribution, which random.expovariate makes. Its parameter is how many arrive per minute, 25 ÷ 60:

import randomimport statistics random.seed(3)gaps = [random.expovariate(25 / 60) for i in range(1000)]print("Average gap:", round(statistics.mean(gaps), 2), "minutes")print("Shortest:", round(min(gaps), 2), "longest:", round(max(gaps), 1))print("Gaps under half a minute:", len([g for g in gaps if g < 0.5]), "of 1000")

Average gap: 2.37 minutesShortest: 0.0 longest: 15.6Gaps under half a minute: 195 of 1000

The gaps average about 2.4 minutes, as they should, but many people come almost together. Those bunches are what build a line, even when the counter could serve everyone if they came evenly.

A simulation of the ticket counter

Now the counter itself. The rules, the model of the counter, are:

• The counter opens at 8:00 and takes new people until 11:00 (180 minutes); everyone in line by then is served.
• People arrive at random, on average per_hour an hour (the exponential gaps).
• Serving one person takes between 1 and 3 minutes, any time in between equally likely (random.uniform(1, 3)), 2 minutes on average. So one counter can serve about 30 people an hour.
• There is one line. The person at its front goes to whichever counter is free first.

Every rule is a simplification of the real counter, and you will come back to that. Here is the model as a function in a module:

"""A simulation of a ticket counter's morning, 8:00 to 11:00.""" import random  def one_morning(per_hour, counters):    """One made-up morning. Gives back each person's wait in line, in minutes."""    free_at = [0] * counters      # when each counter is next free (minutes after 8:00)    waits = []    now = 0    while True:        now = now + random.expovariate(per_hour / 60)   # the next person arrives        if now > 180:                                   # 11:00: no new people            break        k = free_at.index(min(free_at))                 # the counter free first        start = max(now, free_at[k])                    # served when it is free        waits.append(start - now)        free_at[k] = start + random.uniform(1, 3)       # serving takes 1 to 3 minutes    return waits

The heart of it is start = max(now, free_at[k]): a person is served at the moment they arrive, or the moment the counter is free, whichever is later. The wait is the difference. Run one morning with 25 people an hour and one counter:

import randomimport statisticsimport counter random.seed(5)waits = counter.one_morning(25, 1)print(len(waits), "people came")print("Average wait:", round(statistics.mean(waits), 1), "minutes")print("Longest wait:", round(max(waits), 1), "minutes")

73 people cameAverage wait: 5.9 minutesLongest wait: 18.6 minutes

One run is a story; many runs are evidence

Was that a typical morning? You can’t tell from one. Another seed gives another morning, perhaps much better, perhaps much worse. So a simulation is run many times, and the results are looked at together, as data:

import randomimport statisticsimport matplotlib.pyplot as pltimport counter random.seed(5)averages = []for day in range(100):    waits = counter.one_morning(25, 1)    averages.append(statistics.mean(waits)) print("Typical morning (median):", round(statistics.median(averages), 1), "minutes")print("Best morning:", round(min(averages), 1), "worst:", round(max(averages), 1))print("Mornings with an average wait over 10 minutes:", len([a for a in averages if a > 10])) plt.hist(averages, bins=15)plt.xlabel("Average wait in a morning (minutes)")plt.ylabel("Mornings")plt.title("100 simulated mornings, 25 people an hour, one counter")plt.show()

Typical morning (median): 3.4 minutesBest morning: 0.8 worst: 16.8Mornings with an average wait over 10 minutes: 3

The histogram shows how the mornings spread: most are bunched at the low end, a few stretch far to the right. One counter for 25 people an hour gives a typical wait of a few minutes, and an occasional bad morning much longer.

Now the question the hospital asked. Try every number of people with one and two counters:

import randomimport statisticsimport counter random.seed(5)print("per hour | 1 counter | 2 counters   (typical average wait, minutes)")for per_hour in [15, 20, 25, 28, 35, 45]:    row = []    for counters in [1, 2]:        averages = [statistics.mean(counter.one_morning(per_hour, counters)) for day in range(100)]        row.append(statistics.median(averages))    print(f"{per_hour:8} | {row[0]:9.1f} | {row[1]:10.1f}")

per hour | 1 counter | 2 counters   (typical average wait, minutes)      15 |       1.0 |        0.1      20 |       1.7 |        0.1      25 |       3.7 |        0.2      28 |       5.7 |        0.3      35 |      19.6 |        0.5      45 |      48.3 |        1.2

(f"{per_hour:8}" is an f-string with a width: the number is printed in 8 places, so the columns line up; :9.1f is 9 places with one decimal place.)

Read down the first column. The waits don’t grow evenly. From 15 to 28 people an hour they grow slowly, a minute or two at a time. Past 30 an hour, more people come than one counter can serve, and the line grows all morning: a typical wait of about 20 minutes at 35, nearly 50 at 45. With two counters, waits stay around a minute or less all the way to 45. A line grows suddenly as the people arriving get close to what the counter can serve: that is the most useful thing this simulation shows, and it is true of bank counters, petrol pumps, toll gates and the servers of a website.

Is the simulation right?

A simulation can only be as good as its model. Before anyone acts on one, it must be evaluated: what does it answer, how well does it match the real world, and what does it leave out?

1. Check the parameters against real data. Are 25 people an hour, and 1 to 3 minutes each, true? Someone should stand at the counter on a few mornings and write down when each person arrives and how long each takes.
2. Check the results against the real thing. Run the simulation with today’s real settings (one counter, the real number of people) and compare its waits with the waits people really had. If the model says 5 minutes and people really wait 40, something in the model is wrong.
3. List what it leaves out. Our model assumes people arrive evenly through the morning, but many come at 8:00, before the counter opens, so the real line starts long. It assumes nobody gives up and goes home. It assumes every clerk is equally quick, and never takes a break. It ignores patients who are too ill to wait. Each of these could change the answer.
4. Say how sure. Run many mornings, and report the spread (the bad mornings too), not one number.

A chart or a simulation is only the middle of the work. At the end someone must say what it shows, and that is where most mistakes are made: not in the programs, but in the sentences written about their results.

What it answers, and what it doesn’t

Before you believe a chart or a simulation, your own or someone else’s, ask:

• What question does it answer? The districts chart answers “do higher districts get more rain?” for these 77 districts, as they were measured. It doesn’t say why, and it doesn’t answer “will it rain more in my village?”
• Where did the data come from? districts.csv’s rain is not from rain gauges in every district. It comes from a climate data set that estimates the rain everywhere on the map from weather records and a computer model of how mountains steer the clouds, averaged over 1981 to 2010. That is a good estimate, but an estimate, and a district’s rain changes from year to year.
• What does it leave out? A simulation leaves out whatever its model leaves out; a chart, whatever wasn’t collected. A chart of districts’ averages hides the difference between a district’s valleys and its peaks.
• How sure is it? Is it based on 13 rows or 13,000? Does it show the spread, or only an average?
• Could something else explain it? Two things that rise together may both be caused by a third. Grade 7 met this as “going together isn’t causing”; the next section shows how far it can go.

A third variable can turn a result round

Here is a puzzle, with made-up numbers. A big city hospital and a district hospital each treated 500 patients last year. At the district hospital 93% recovered; at the city hospital only 84%. Which hospital is better?

The district hospital, surely. But split the patients by how ill they were when they came (Figure 4.5):

In both groups the city hospital does better, yet over all its patients it does worse. How? The city hospital takes most of the seriously ill patients, the ones sent on from other hospitals, and they recover less often wherever they are treated. Its overall rate is low because of who it treats, not how well. A result that holds in every group but turns round when the groups are added together is called Simpson’s paradox. The cure is the multivariate thinking of this unit: before comparing, ask what third variable (here, how ill the patients were) differs between the groups, and split the data by it.

Writing a conclusion

A good conclusion has four parts: the claim, the evidence for it, how sure you are, and the limits. It uses careful words: “in these data”, “tend to”, “about”, “we found no evidence that…”, instead of “always”, “proves” or “everyone”.