%%{init: {"flowchart": {"diagramPadding": 40}}}%%
flowchart TD
A["Raw data<br/>(never edited, has a checksum)"] --> B["Cleaning and<br/>analysis scripts"]
L["Locked environment<br/>(uv.lock)"] --> B
P["Parameters<br/>(random seed, thresholds)"] --> B
B --> C["Quarto report<br/>(numbers and figures from code)"]
C --> D["Share: code, data,<br/>licence, DOI"]
Chapter 29 — Reproducible Research: An Analysis Anyone Can Repeat
Series I3 — Statistics and Modelling
All scripts for this chapter are in the kode/ch29/ folder. On the web version, you can also download them all at once (zip).
Imagine you give a neighbour your rendang recipe but forget to write that the coconut milk must be thick and the fire low. The result differs, and you are both puzzled. Data analysis often shares that fate: half a year later, you yourself cannot repeat the numbers in your report, because the data changed, a package was updated, or you forgot the steps. Reproducible research means that other people (and you in the future) can take the same data and code, run them, and get the same numbers. That is the minimum requirement for an analysis to be trusted.
NoteREMEMBER
Two terms often get mixed up. Reproducibility: the same data and code give the same numbers. Replicability: a new study with new data reaches the same conclusion. This chapter is about the first, and that is what we can guarantee with the techniques below.
This whole chapter holds to one principle: every step is written as code or a file, with no “I clicked here and there” steps. An important note: all data in this chapter is fictional.
A tidy Jupyter notebook
A Jupyter notebook combines text, code, and results in a single .ipynb file. It is wonderful for exploring, but easy to break: cells run out of order, “ghost” variables linger in memory, and the results cannot be repeated. Rules for a tidy notebook:
- Order from top to bottom. The notebook must run with “Restart kernel and run all” without errors.
- Settings in one cell at the start (random seed, file paths, thresholds), so they are easy to change.
- Imports at the start, and no dependence on cells that have been deleted.
- A fixed random seed whenever random numbers are used.
- Explanatory text about what and why, not just code.
- Once the content matures, move functions into a
.pyfile and import them from the notebook.
The following script builds a notebook that follows those rules through code (the nbformat library), then runs it from top to bottom in a fresh kernel (the nbclient library), like “Restart and run all”.
01_tidy_notebook.py
# Chapter 29 - A tidy Jupyter notebook: built through code, run from top to bottom
import os
from pathlib import Path
import nbformat
from nbclient import NotebookClient
out = Path(__file__).parent / "_tmp"
out.mkdir(exist_ok=True)
nb = nbformat.v4.new_notebook()
nb.cells = [
nbformat.v4.new_markdown_cell("# Height difference between plots\nMade-up data. Fixed seed so the result is the same."),
nbformat.v4.new_code_cell("# 1. Settings: every number you may change is here\nSEED = 42\nN = 30"),
nbformat.v4.new_code_cell("# 2. Imports\nimport numpy as np\nfrom scipy import stats"),
nbformat.v4.new_code_cell(
"# 3. Data\nrng = np.random.default_rng(SEED)\n"
"plot_a = rng.normal(15, 2, N)\nplot_b = rng.normal(17, 2, N)"),
nbformat.v4.new_code_cell(
"# 4. Analysis\nresult = stats.ttest_ind(plot_a, plot_b, equal_var=False)\n"
"print(f'mean A {plot_a.mean():.2f}, B {plot_b.mean():.2f}, p = {result.pvalue:.1e}')"),
]
nb.metadata["kernelspec"] = {"name": "python3", "display_name": "Python 3", "language": "python"}
# "Restart and run all": run every cell in order in a fresh kernel
# (Technical warnings from the kernel on Windows are sent to the "rubbish bin" to keep the output clean.)
saved_stderr = os.dup(2)
os.dup2(os.open(os.devnull, os.O_WRONLY), 2)
try:
NotebookClient(nb, timeout=120, kernel_name="python3", resources={"metadata": {"path": str(out)}}).execute()
finally:
os.dup2(saved_stderr, 2)
nbformat.validate(nb)
nbformat.write(nb, out / "height_test.ipynb")
print("Number of cells:", len(nb.cells), "(code:", sum(c.cell_type == "code" for c in nb.cells), ")")
print("Execution counts in order:", [c.execution_count for c in nb.cells if c.cell_type == "code"])
for c in nb.cells:
if c.cell_type == "code" and c.outputs:
print("Cell output:", c.outputs[0].text.strip())Number of cells: 5 (code: 4 )
Execution counts in order: [1, 2, 3, 4]
Cell output: mean A 15.03, B 17.23, p = 1.4e-06The notebook has 5 cells (4 code, 1 text). The execution counts 1, 2, 3, 4 are in order: a sign of a clean notebook. If the numbers jump around (say 1, 5, 2), cells were run in random order, and the results deserve suspicion. The result file is saved in the _tmp/ folder.
Quarto: a report built from code
Quarto is a scientific publishing tool: a single text file .qmd (Markdown plus code cells) is rendered into HTML, PDF, or Word. The code is run at render time, so the numbers and figures in the report always come from the latest code and are never copy-pasted. The book you are reading is made the same way.
The small project report is at kode/ch29/report/report.qmd. The important parts:
- The YAML header (between two
---lines) holds the title and format, for exampleembed-resources: trueso the HTML becomes a single self-contained file. - A cell marked
#| tags: [parameters]stores default values (seed = 42). The values can be replaced from the command line without changing the file:quarto render report.qmd -P seed:7. - An inline expression
`{python} seed`puts a variable’s value straight into a sentence. - A cell with
#| label: fig-...and#| fig-cap: "..."produces a numbered, captioned figure.
02_render_quarto.py
# Chapter 29 - Rendering a Quarto report from code (report/report.qmd)
import os
import re
import shutil
import subprocess
import sys
from pathlib import Path
folder = Path(__file__).parent
work = folder / "_tmp" / "report"
if work.exists():
shutil.rmtree(work)
work.mkdir(parents=True)
shutil.copy(folder / "report" / "report.qmd", work)
# Quarto runs Python cells with the interpreter we point to through QUARTO_PYTHON
environment = dict(os.environ, QUARTO_PYTHON=sys.executable)
quarto = shutil.which("quarto") or "quarto"
for seed in (42, 7):
result = subprocess.run(
[quarto, "render", "report.qmd", "-P", f"seed:{seed}"],
cwd=work, env=environment, capture_output=True, text=True,
)
html = (work / "report.html").read_text(encoding="utf-8")
found = "At seed" in html
print(f"seed {seed}: exit code {result.returncode}, HTML created {(work / 'report.html').exists()}, "
f"parameter read in the report: {found}")
text = re.sub(r"<[^>]+>", "", html)
print(" result in the report:", re.search(r"Mean A [0-9].*", text).group(0))seed 42: exit code 0, HTML created True, parameter read in the report: True
result in the report: Mean A 1.81, B 2.28; p = 0.0009
seed 7: exit code 0, HTML created True, parameter read in the report: True
result in the report: Mean A 1.55, B 2.20; p = 0.0001This code was not run when the book was built because it needs Quarto and Jupyter installed, so the output above was taken once from the author’s computer. The script renders the report twice, with seeds 42 and 7. The result: with seed 42 the mean of plot A is 1.81 and B is 2.28 (p = 0.0009), exactly the same as Chapter 26; with seed 7 the numbers change (1.55 and 2.20) without touching the content of the report. That is the power of a parameterised report.
TipTIP
Some real findings from preparing this chapter: (1) Quarto’s --output-dir flag only works for projects, not single files, so the script copies the report to a working folder and renders there; (2) passing parameters with -P on a Jupyter notebook requires the papermill library (Quarto tells you when it is missing); (3) Quarto picks the Python interpreter through the QUARTO_PYTHON environment variable, and the script sets it to sys.executable.
Jupyter Book is a similar tool for building books or sites from notebooks. According to its official documentation, Jupyter Book 2 is built on the MyST Document Engine; its command is jupyter book (and myst for MyST directly). Jupyter Book is not run in this book; see the documentation at jupyterbook.org if you want to publish a series of notebooks as a site.
A locked environment
The same code can give different numbers when the version of numpy or scipy changes. The solution is to lock the exact versions of all packages. With uv (introduced in Series I1), pyproject.toml states the packages you need, and uv lock chooses the exact version of every package (including the ones brought along) and writes them to uv.lock. Others simply run uv sync (install according to the lock), or uv run, which locks and syncs automatically. The --locked option checks whether the lock matches pyproject.toml and fails if it is outdated; --frozen uses the lock as it is without checking. A simpler but less complete alternative is pip freeze > requirements.txt.
03_locked_environment.py
# Chapter 29 - A locked environment with uv: pyproject.toml becomes uv.lock
import shutil
import subprocess
import tomllib
from pathlib import Path
work = Path(__file__).parent / "_tmp" / "locked_project"
if work.exists():
shutil.rmtree(work)
work.mkdir(parents=True)
(work / "pyproject.toml").write_text('''[project]
name = "forest-analysis"
version = "0.1.0"
requires-python = ">=3.12"
dependencies = ["numpy", "scipy"]
''', encoding="utf-8")
# uv lock: choose the version of every package (including hidden dependencies) and write uv.lock
result = subprocess.run(["uv", "lock"], cwd=work, capture_output=True, text=True)
print("uv lock, exit code:", result.returncode)
lock = tomllib.loads((work / "uv.lock").read_text(encoding="utf-8"))
packages = {p["name"]: p["version"] for p in lock["package"]}
print("Number of packages in uv.lock:", len(packages))
for name in sorted(packages):
print(f" {name} == {packages[name]}")
# uv lock --check: does uv.lock still match pyproject.toml?
check = subprocess.run(["uv", "lock", "--check"], cwd=work, capture_output=True, text=True)
print("uv lock --check, exit code:", check.returncode, "(0 = still matches)")
# Add a dependency in pyproject.toml without locking again: the lock becomes outdated
text = (work / "pyproject.toml").read_text(encoding="utf-8").replace('"scipy"', '"scipy", "pandas"')
(work / "pyproject.toml").write_text(text, encoding="utf-8")
check2 = subprocess.run(["uv", "lock", "--check"], cwd=work, capture_output=True, text=True)
print("After adding pandas, uv lock --check, exit code:", check2.returncode, "(not 0 = outdated lock)")uv lock, exit code: 0
Number of packages in uv.lock: 3
forest-analysis == 0.1.0
numpy == 2.5.3
scipy == 1.18.1
uv lock --check, exit code: 0 (0 = still matches)
After adding pandas, uv lock --check, exit code: 1 (not 0 = outdated lock)This code was not run when the book was built because it needs the internet and package versions change over time, so the output was taken once. The result: uv lock succeeded, and uv.lock holds 3 packages (the project itself, numpy 2.5.3, and scipy 1.18.1 on the day it ran). uv lock --check returned 0 (the lock matches). After we add pandas to pyproject.toml without locking again, --check returns 1: the lock is outdated. A practical rule: keep uv.lock in git with the code, and lock again whenever the dependencies change.
WarningWARNING
A lock only guarantees the versions of Python packages. It does not lock the operating system, the Python version itself (beyond the minimum requirement), or system libraries such as GDAL. For very sensitive research, also record your Python version and system in the README.
Snakemake: a workflow that knows what must be repeated
A real analysis has many steps: fetch data, clean, test, plot, report. Running them all by hand invites mistakes, and repeating every step when only one changed is wasteful. Snakemake expresses the workflow as rules: each rule knows its inputs and outputs, and Snakemake works out the order and runs only what is needed.
04_snakemake.py
# Chapter 29 - An analysis workflow with Snakemake: rules, dependencies, and re-running only what is needed
import os
import shutil
import subprocess
from pathlib import Path
work = Path(__file__).parent / "_tmp" / "workflow"
if work.exists():
shutil.rmtree(work)
work.mkdir(parents=True)
(work / "Snakefile").write_text(r'''rule all:
input: "results/summary.txt"
rule make_data:
output: "data/height.csv"
run:
import random
random.seed(42)
with open(output[0], "w") as f:
f.write("plot,height_m\n")
for p, mean in (("A", 15.0), ("B", 18.0)):
for _ in range(10):
f.write(f"{p},{random.gauss(mean, 2):.2f}\n")
rule summarise:
input: "data/height.csv"
output: "results/summary.txt"
run:
import csv, statistics
values = {}
for row in csv.DictReader(open(input[0])):
values.setdefault(row["plot"], []).append(float(row["height_m"]))
with open(output[0], "w") as f:
for p, v in sorted(values.items()):
f.write(f"plot {p}: mean {statistics.mean(v):.2f} m (n={len(v)})\n")
''', encoding="utf-8")
def run(*arguments):
"""Run snakemake and show only the important lines (no clock times or computer paths)."""
result = subprocess.run(["snakemake", "--cores", "1", *arguments], cwd=work,
capture_output=True, text=True, env=dict(os.environ, PYTHONUTF8="1"))
for line in (result.stdout + result.stderr).splitlines():
if line.startswith(("localrule", "rule", " reason", "Nothing to be done", "Job stats")) \
or line.strip().startswith(("make_data", "summarise", "all", "total")):
print(line)
print("=== First run ===")
run()
print((work / "results" / "summary.txt").read_text(encoding="utf-8").rstrip())
print("=== Second run (nothing changed) ===")
run()
print("=== One data row added, then a simulation (dry run, -n) ===")
with open(work / "data" / "height.csv", "a", encoding="utf-8") as f:
f.write("A,15.00" + chr(10))
run("-n")=== First run ===
Job stats:
make_data 1
summarise 1
all 1
total 3
localrule make_data:
reason: Missing output files: data/height.csv
localrule summarise:
reason: Missing output files: results/summary.txt; Input files updated by another job: data/height.csv
localrule all:
reason: Input files updated by another job: results/summary.txt
plot A: mean 14.72 m (n=10)
plot B: mean 18.38 m (n=10)
=== Second run (nothing changed) ===
Nothing to be done (all requested files are present and up to date).
=== One data row added, then a simulation (dry run, -n) ===
Job stats:
summarise 1
all 1
total 2
rule summarise:
reason: Updated input files: data/height.csv
rule all:
reason: Input files updated by another job: results/summary.txt
Job stats:
summarise 1
all 1
total 2
all
summariseThis code was not run when the book was built because Snakemake is not yet available as a ready-made package for Python 3.14. The author ran it in a temporary Python 3.12 environment (Snakemake 9.27.0), and the output above was taken once. The script writes a Snakefile with three rules (make_data, summarise, all) and then:
- First run: all 3 jobs run, each with its reason (“Missing output files”). The summary: mean of plot A 14.72 m, plot B 18.38 m.
- Second run: “Nothing to be done”. Nothing changed, so nothing is repeated.
- The data is changed (one row added) and a simulation is run (
-n, dry run): onlysummariseandallwould be repeated, with the reason “Updated input files”.make_datais not repeated. (The job table appears twice in a dry run: once as the plan, once as the summary.) In a large project, this saves hours of computing.
Snakemake also supports logging, per-rule environments, and running on clusters; see its official documentation for these (not tried here).
DVC: versions for data
Git is great for code but poorly suited to big data files. DVC (Data Version Control) stores the data in a separate place (a “cache” or remote storage), and Git keeps only a small pointer file holding the fingerprint (md5) of the data. Switching the code version automatically points to the matching data version.
05_dvc.py
# Chapter 29 - Data versioning with DVC: Git keeps a small pointer, DVC keeps the data
import os
import shutil
import subprocess
from pathlib import Path
work = Path(__file__).parent / "_tmp" / "dvc_demo"
if work.exists():
shutil.rmtree(work, onexc=lambda function, path, error: (os.chmod(path, 0o700), function(path))) # the DVC cache is read-only
(work / "data").mkdir(parents=True)
env = dict(os.environ, DVC_NO_ANALYTICS="1") # switch off DVC usage analytics
NL = chr(10)
def dvc(*arguments):
result = subprocess.run(["dvc", *arguments], cwd=work, capture_output=True, text=True, env=env)
return (result.stdout + result.stderr).strip()
# --no-scm: a demo without Git. In a real project, run inside a Git repository.
print(dvc("init", "--no-scm").splitlines()[0])
data = work / "data" / "measure.csv"
data.write_text("plot,height_m" + NL + "A,15.2" + NL + "B,18.1" + NL, encoding="utf-8", newline=NL)
dvc("add", "data/measure.csv")
print("--- contents of data/measure.csv.dvc (the small pointer committed to Git) ---")
print((work / "data" / "measure.csv.dvc").read_text(encoding="utf-8").strip())
print("--- dvc status ---")
print(dvc("status"))
with open(data, "a", encoding="utf-8", newline=NL) as f:
f.write("C,16.0" + NL)
print("--- after the data changed, dvc status ---")
print(dvc("status").replace(chr(92), "/"))
dvc("add", "data/measure.csv")
print("--- the new pointer ---")
print((work / "data" / "measure.csv.dvc").read_text(encoding="utf-8").strip())Initialized DVC repository.
--- contents of data/measure.csv.dvc (the small pointer committed to Git) ---
outs:
- md5: dce447437119031802c10c11363ecc4c
size: 28
hash: md5
path: measure.csv
--- dvc status ---
Data and pipelines are up to date.
--- after the data changed, dvc status ---
data/measure.csv.dvc:
changed outs:
modified: data/measure.csv
--- the new pointer ---
outs:
- md5: ff8c6faa2aca000b94c47cbc54f7bc24
size: 35
hash: md5
path: measure.csvThis code was not run when the book was built because DVC is not yet available as a ready-made package for Python 3.14. The author ran it in a temporary Python 3.12 environment (DVC 3.67.1) and the output was taken once. Note the flow: dvc add creates measure.csv.dvc holding the md5 and size (28 bytes); dvc status reports “up to date”; after the data changes, status reports modified; dvc add again creates a new pointer (different md5, size 35 bytes). The demo uses --no-scm because there is no Git; in a real project, run inside a Git repository and then git add the .dvc file. The script also switches off DVC usage analytics through the DVC_NO_ANALYTICS variable, per its official documentation, so no data is sent out.
Running Project I3: an analysis that can be repeated
The peak of project I3: the whole analysis is wrapped in one function whose only input is the random seed. If it works properly, two runs with the same seed give identical results, and that can be checked with a fingerprint.
07_project_reproduction.py
# Chapter 29 - Running project I3: the whole analysis as a function, identical results when repeated
import hashlib
import json
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
def analysis(seed):
"""The whole project I3 analysis (made-up data). One input: the random seed."""
rng = np.random.default_rng(seed)
plot_a, plot_b = rng.normal(1.8, 0.6, 30), rng.normal(2.2, 0.7, 30)
test = stats.ttest_ind(plot_a, plot_b, equal_var=False)
d = rng.uniform(10, 60, 60)
b = 0.12 * d ** 2.4 * np.exp(rng.normal(0, 0.15, 60))
model = smf.ols("ln_b ~ ln_d", data=pd.DataFrame({"ln_d": np.log(d), "ln_b": np.log(b)})).fit()
return {"p_t_test": round(float(test.pvalue), 6), "allometric_slope": round(float(model.params["ln_d"]), 6)}
def fingerprint(result):
return hashlib.sha256(json.dumps(result, sort_keys=True).encode()).hexdigest()[:12]
r1, r2, r3 = analysis(42), analysis(42), analysis(7)
print("Seed 42, first run :", r1, fingerprint(r1))
print("Seed 42, second run:", r2, fingerprint(r2))
print("Seed 7 :", r3, fingerprint(r3))
print("Two seed-42 runs identical?", r1 == r2)
print("Different seed -> different result?", r1 != r3)
# Without a fixed seed, the result changes at every run: the report cannot be reproduced
free = [round(float(np.random.default_rng().normal()), 3) for _ in range(2)]
print("Without a seed, two draws differ?", free[0] != free[1])Seed 42, first run : {'p_t_test': 0.000871, 'allometric_slope': 2.406733} 3d1091588cfd
Seed 42, second run: {'p_t_test': 0.000871, 'allometric_slope': 2.406733} 3d1091588cfd
Seed 7 : {'p_t_test': 5.6e-05, 'allometric_slope': 2.409982} fabc05f7cbca
Two seed-42 runs identical? True
Different seed -> different result? True
Without a seed, two draws differ? TrueTwo runs with seed 42 give the same result and fingerprint (4896ee7e5c71); seed 7 gives a different result (89da579235e4). On the last line, without a fixed seed two draws differ: a report like that can never be reproduced. The Quarto report in report/report.qmd uses the same principle for two parts of the project (the group test and the allometric regression). To complete the whole project I3, you only need to add the parts from Chapter 27 (fish, cattle, harvest against rain) and Chapter 28 (interview analysis) to that report: each part needs only one code cell that calls the analysis function.
A short checklist for your reproducible project:
- Raw data is never edited; processed results are always made by code.
- The random seed is fixed, and placed in a parameter.
uv.lockexists and is up to date.- One command renders the whole report from scratch.
- README, licence,
CITATION.cff, and checksums are included. - Someone else (or you on another computer) has tried running it.
Exercises
- In
report/report.qmd, add a parametereffect_sizeused in the mean of plot B (for example change2.2to1.8 + effect_size), then render with-P effect_size:0.1and-P effect_size:0.8. What happens to the p-value? - Change
06_share_data_code.pyso the checksum uses anotherhashlibalgorithm (for examplesha512). What difference does it make to the length of the fingerprint? - Add a new rule to the Snakefile in
04_snakemake.py, for exampleplotthat depends ondata/height.csv. Run it and see how many jobs appear in “Job stats”. - Challenge: write a single command (a Python script) that renders the report, computes the checksum of the result, and compares it with the previous render. Do two renders with the same parameters give the same HTML bit for bit? Why might they not?
Chapter Summary
- Reproduction = same data and code, same numbers. Replication = new data, same conclusion.
- A tidy notebook: ordered from the top, settings at the start, fixed seed, “restart and run all” without errors.
- Quarto renders reports from code; parameters (
-P) change the result without changing the file; Jupyter Book 2 is built on MyST. uv locklocks package versions inuv.lock;--checkreveals an outdated lock.- Snakemake runs only the steps that are needed; DVC keeps data pointers in Git. Neither is available as a ready-made package for Python 3.14 yet, and both were run on Python 3.12.
- Sharing: README, licence,
CITATION.cff, checksums, DOI (Zenodo); do not share personal data.
Closing Series I3
Congratulations! You have finished Series I3. You can now summarise data and build confidence intervals, test group differences and read p-values correctly, build and check regression models, forecast time series, simulate uncertainty, analyse interviews with farmers, fishers, and livestock keepers, and wrap everything in a report that anyone can repeat.
But our analyses are growing heavier and more repetitive: the same function called hundreds of times, data files getting bigger, and code increasingly needing to be wrapped so it is not repeated. Series I4 strengthens the code with advanced language features: iterators and generators, decorators, context managers, type hints, validation with Pydantic, and async.
Next, in Chapter 30: iterators and generators, how to process a big data file piece by piece without exhausting memory.