L1 demo: an analysis that does not reproduce

L1 demo: an analysis that does not reproduce#

This notebook loads a year of roadside air-quality sensor data, plots a reference measurement against the corresponding sensor channel, fits a naive calibration, and integrates the concentration curve to get a cumulative exposure.

All of that is ordinary work. The notebook is nonetheless broken in three independent ways, and you will meet them in order as you run it from the top.

Your job is to name each one before I do, and to say not just what the fix is but what practice would have stopped it reaching a colleague. None of the three are exotic. That is the point: reproducibility fails for boring reasons, and boring reasons are the ones that survive review.

Data: UCI Air Quality Data Set, De Vito et al. Hourly readings from a multisensor device deployed on an Italian roadside, March 2004 onward.

Run this first on Colab#

Colab starts from its own preinstalled environment rather than this course’s uv environment, so run the cell below before anything else. It installs what this notebook needs and Colab does not already have. Outside Colab it does nothing, so you can run it or skip it.

# Run this first on Colab. Anywhere else this cell does nothing.
#
# Only genuinely missing packages are installed, so Colab's own versions of
# everything it already ships are left alone.
import importlib.util
import subprocess
import sys

REQUIREMENTS = {
    "matplotlib": "matplotlib",
    "numpy": "numpy",
    "pandas": "pandas",
    "sklearn": "scikit-learn",
}


def _missing(module):
    try:
        return importlib.util.find_spec(module) is None
    except ModuleNotFoundError:  # the parent package is absent
        return True


if "google.colab" in sys.modules:
    need = sorted({pip for mod, pip in REQUIREMENTS.items() if _missing(mod)})
    if need:
        print("installing:", " ".join(need))
        subprocess.run([sys.executable, "-m", "pip", "install", "-q", *need], check=True)
    print("Colab setup done." if need else "Colab: nothing to install.")

1. Environment#

These imports pull in NumPy, pandas, and scikit-learn, but nothing here records which version of any of them this analysis was written against, and there is no lockfile beside the notebook.

Print the versions and hold onto them. They matter in section 5.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score

print('numpy ', np.__version__)
print('pandas', pd.__version__)

2. Load the data#

The raw CSV is not in this repository. Datasets do not belong in git: they bloat the history, they are usually someone else’s to license, and a repository is not a distribution channel. So the next cell fetches it from UCI on first run and reuses the local copy afterwards, which is what you should do in your own projects too.

The export has the usual instrument-file quirks, none of them announced: semicolon separated, comma as the decimal mark, two trailing empty columns, and missing values coded as -200 rather than left blank.

import io
import urllib.request
import zipfile
from pathlib import Path

DATA_PATH = Path('/Users/jkitchin/Dropbox/classes/f26-systems-toolchains/data/AirQualityUCI.csv')
URL = 'https://archive.ics.uci.edu/static/public/360/air+quality.zip'

if not DATA_PATH.exists():
    print(f'fetching {URL}')
    with urllib.request.urlopen(URL) as response:
        payload = response.read()
    with zipfile.ZipFile(io.BytesIO(payload)) as archive:
        DATA_PATH.write_bytes(archive.read('AirQualityUCI.csv'))

print(f'using {DATA_PATH}')

If that cell just failed#

Read the error. It has nothing to do with statistics, sensors, or modelling. It is a FileNotFoundError, and it happened while trying to write the downloaded file.

The download worked. The place it was told to put the file does not exist on your machine, and never will.

Problem one. Fix it, rerun, and the fetch will succeed. Then keep going.

df = (
    pd.read_csv(DATA_PATH, sep=';', decimal=',')
    .dropna(axis=1, how='all')
    .dropna(how='all')
)

df['ts'] = pd.to_datetime(
    df['Date'] + ' ' + df['Time'].str.replace('.', ':', regex=False),
    format='%d/%m/%Y %H:%M:%S',
)
df = df.replace(-200, np.nan)

print(f'{len(df)} rows, {df.ts.min().date()} to {df.ts.max().date()}')
df.head(3)

3. A first look#

CO(GT) is a reference-grade carbon monoxide measurement. PT08.S1(CO) is the cheap metal-oxide sensor that is supposed to track it. If the sensor is any good, these should be related.

d = df.dropna(subset=['CO(GT)', 'PT08.S1(CO)'])

fig, ax = plt.subplots(figsize=(6, 4.2))
ax.scatter(d['PT08.S1(CO)'], d['CO(GT)'], s=4, alpha=0.25)
ax.set_xlabel('PT08.S1(CO) sensor response')
ax.set_ylabel('CO(GT) reference, mg/m$^3$')
ax.set_title('Sensor response vs reference measurement')
plt.show()

4. A naive calibration#

Fit a linear model from the sensor channel to the reference, hold out a quarter of the data, and score it.

X = d[['PT08.S1(CO)']]
y = d['CO(GT)']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)

model = LinearRegression().fit(X_train, y_train)
score = r2_score(y_test, model.predict(X_test))

print(f'R2 on held-out data: {score:.4f}')
print(f'slope {model.coef_[0]:.5f}, intercept {model.intercept_:.5f}')

Run that cell again#

Do not change anything. Run it a second time, and a third. Write the numbers down.

Then decide which of them you would be willing to put in a report, and how you would answer a reviewer who asked you to reproduce it six months from now.

Problem two.

5. Cumulative exposure#

A calibration is not the quantity a health question actually needs. What matters is cumulative dose: the area under the concentration curve over the deployment. That is a numerical integration over the hourly series.

We do it for two pollutants. Run both cells.

hourly_co = d.set_index('ts')['CO(GT)'].resample('h').mean().dropna()
hours_co = (hourly_co.index - hourly_co.index[0]).total_seconds() / 3600.0

exposure_co = np.trapz(hourly_co.values, hours_co)
print(f'Integrated CO exposure:  {exposure_co:,.0f} mg-h/m3')

nox = df.dropna(subset=['NOx(GT)'])
hourly_nox = nox.set_index('ts')['NOx(GT)'].resample('h').mean().dropna()
hours_nox = (hourly_nox.index - hourly_nox.index[0]).total_seconds() / 3600.0

exposure_nox = np.trapezoid(hourly_nox.values, hours_nox)
print(f'Integrated NOx exposure: {exposure_nox:,.0f} ppb-h')

Exactly one of those two cells failed#

Not both, and not neither. Which one depends on the NumPy version you printed in section 1.

Compare with the person next to you. If you got different failures, neither of you did anything wrong, and that is the entire lesson. Both cells are ordinary numerical integration. Both were correct when someone wrote them.

Problem three. Note that nothing in this repository told you which version this analysis expected, so nothing told you which cell was the broken one.

Your turn#

Fill this in. The last column is the one that matters.

#

Symptom

Root cause

Fix

Practice that prevents it

1

2

3

Two questions to take away:

  • Problem two does not raise an error. Nothing in the output marks it as wrong. How many analyses have you run that had this property?

  • Problem three was introduced by a library doing something entirely reasonable, a deprecated name being removed in a major release. Whose responsibility was it to notice?

We rebuild this properly in L2: a uv-managed project with a locked environment, a pinned interpreter, relative paths, seeded splits, and the analysis moved out of the notebook into a module you can test and rerun. That is also assignment A1.