Statistical software and reproducibility — Regression & Experimental Design

What you'll learn

Rerun an analysis in R or Python
Compare the main statistical tools and know that their defaults can differ
Apply reproducibility habits: read-only raw data, one script, seeds, versions and version control
Cross-check an important result in a second tool

Compare R and Python for statistics and rerun an analysis reproducibly in a second tool.

Why this matters

A statistical result is only as trustworthy as your ability to produce it again. Reviewers, auditors and future colleagues will ask for it, and so will you, six months later, when new data arrive or a number in the report looks odd. Reproducibility is also how errors get caught: a famous 2013 re-analysis of an influential economics paper found a spreadsheet error simply because the original authors shared their file. This lesson covers the main statistical software tools and the habits that make an analysis rerunnable.

The main tools

Tool	Strengths	Typical users
R	Widest range of statistical methods; `lme4`, `emmeans`, `ggplot2`, the tidyverse; R Markdown / Quarto	Statisticians, biostatistics, academic research
Python	`pandas`, `statsmodels`, `scipy`, `scikit-learn`, `matplotlib`; one language from data cleaning to production	Data science, ML, engineering, OR
SAS	Long history in pharma and government; validated procedures	Clinical trials, regulated industries
JMP, Minitab, Design-Expert	Point-and-click DOE and quality tools	Process engineers, Six Sigma, labs
Stata, SPSS	Survey and social-science analyses	Economics, public health, social sciences

You do not need all of them. Be fluent in one of R or Python, able to read the other, and comfortable reading output from the point-and-click tools your collaborators use. Most methods in this course exist in both:

Task	R	Python
Linear model	`lm(y ~ x1 + x2)`	`smf.ols("y ~ x1 + x2", df).fit()`
GLM	`glm(y ~ x, family = poisson)`	`smf.glm("y ~ x", df, family=sm.families.Poisson()).fit()`
ANOVA table	`anova(model)` / `car::Anova`	`anova_lm(model, typ=2)`
Mixed model	`lme4::lmer(y ~ x + (1 \| g))`	`smf.mixedlm("y ~ x", df, groups=df["g"]).fit()`
Tukey HSD	`TukeyHSD(aov(...))`	`pairwise_tukeyhsd(y, group)`

Defaults differ. R's anova() gives sequential (Type I) sums of squares; statsmodels' anova_lm(typ=2) gives Type II. For unbalanced data these can disagree. Contrast coding, optimisers and degrees-of-freedom methods for mixed models also differ. When two tools give different answers, it is usually a default, and you should be able to find which one.

Habits of a reproducible analysis

Raw data is read-only. Never edit the original file. All cleaning happens in code, so every change is recorded.
One script (or notebook) runs everything, from raw data to every number and figure in the report. If a number in the report was typed by hand, it can drift from the analysis.
Fix random seeds for randomisation, simulation, bootstrapping and train/test splits, and record them.
Record versions. In Python, pip freeze > requirements.txt or a pyproject.toml; in R, renv::snapshot(). Note the language version too.
Use version control (Git) for code and the report. Commit messages explain why things changed.
Literate documents (Quarto, R Markdown, Jupyter) mix code, output and text so the report regenerates from data. Restart and run all before sharing a notebook; out-of-order cells are a common source of irreproducible results.
Organise the project predictably:

project/
  data/raw/          original files, never edited
  data/processed/    written by code only
  src/ or R/         analysis code
  reports/           Quarto / notebook that builds the report
  requirements.txt   or renv.lock
  README.md          how to rerun, in three lines

Cross-checking in a second tool

A strong check on an important result is to rerun it independently: a colleague reproduces it from the methods section alone, or you reproduce it in a second tool. For the catalyst ANOVA from module 2:

Python

from scipy import stats
print(stats.f_oneway([82, 85, 88, 84, 86], [90, 92, 89, 93, 91], [84, 83, 86, 85, 82]))
# F_onewayResult(statistic=21.5, pvalue=0.000108...)

d <- data.frame(y = c(82,85,88,84,86, 90,92,89,93,91, 84,83,86,85,82),
                catalyst = rep(c("A","B","C"), each = 5))
summary(aov(y ~ catalyst, data = d))   # F value 21.5

Free ways to try both without installing anything: Google Colab for Python and webR or Posit Cloud (free tier) for R. Free tiers and their limits change; check the provider's page before relying on one.

Common mistakes

Copying numbers by hand from output into a report.
"It works on my machine": no record of package versions.
Unseeded randomness, so the bootstrap interval changes every run.
Hidden state in notebooks (cells run out of order, variables defined in deleted cells).
Overwriting raw data while cleaning.
Trusting one tool's default without knowing what it is.

Quick recap

Know R or Python well and read the other; know how their defaults differ.
Raw data read-only, one script end to end, seeds, versions, Git, literate reports.
Cross-check important results in a second tool or by a colleague.

Try it

In the Do tab you will rerun an analysis from this course in a second tool (R or Python), record the environment and confirm the numbers match.

Key terms (5)

Reproducibility: The ability to produce the same result again from the same data and code.
Raw data is read-only: The habit of never editing original data files, so every cleaning step is recorded in code.
Literate documents: Files such as Quarto, R Markdown or Jupyter notebooks that mix code, output and text so a report regenerates from data.
Type II: The kind of sums of squares statsmodels' anova_lm(typ=2) gives, which can differ from R's sequential Type I for unbalanced data.
version control: Tracking changes to code and reports with a tool such as Git, with messages explaining why things changed.