Week 9: Checking Your Work

Reproducibility, reporting, and knowing what AI can’t do for you

Introduction

Your projects are nearly finished. Before you write up, check the work with fresh eyes — because a reviewer will, and so will whoever tries to reproduce it.

Today has four exercises. Two are about your own code being trustworthy. Two are about catching an error in someone else’s, where “someone else” is a language model that will be confident, fluent, and wrong.

Goals:

  • Make a broken script run from a cold start, and know why it broke
  • Report a result in a form a policy reader can use
  • Catch a confident wrong answer produced by plausible-looking code
  • Ask “how plausible was this before we tested?” one last time

The code boxes start empty, except for comments to guide you through the code that’s expected.

Work one comment at a time. If you get stuck:

  1. Hint 1 restates what you are trying to do.
  2. Hint 2 sketches the shape the code might take.
  3. Hint 3 gives you the code with the key pieces blanked out.
  4. Solution gives you the lot.

Take a real attempt before opening a hint: you will learn by doing, not by reading. The demonstrators are here to help.

Three things worth knowing:

  • The boxes share one R session, exactly like a script. Anything you create in Exercise 1 is still there in Exercise 2.
  • Variable names are up to you. The last thing your block prints is what gets checked, so end each block with the answer.
  • Exercises have extensions underneath them. If you finish early, give them a stab. These are designed to equip you to tackle real coding problems you might face, using the documentation. We can help – but check the manuals first!
TipStuck?

Ask on the class board, in Code Q&A. Someone else is likely stuck on the same thing, so asking in public helps them too. Demonstrators check it most weekdays, but they give classmates a chance to answer first.

The board is private to this class, so sign in to GitHub first. If the link shows “404 – page not found”, you’re not signed in, or you haven’t accepted the Classroom 50 invitation yet.

Exercise 1: Make it run from cold

Below is a script that worked perfectly on the machine it was written on, and fails for everyone else. It has four faults, and each is one of the four commonest reasons a shared analysis will not run.

Fix all four so the block runs start to finish and ends by printing a reproducible bootstrap confidence interval.

Use set.seed(2847) so everyone in the room gets the same answer — which is itself the point of the fourth fault.

NoteHint 1

Run it and read the first error. Fix that, run again, read the next. Do not try to spot all four by eye — working through errors one at a time is the actual skill, and it is faster.

The four faults, in the order R will meet them:

  1. A path that exists on one computer in the world.
  2. Two variables defined in a comment, which is to say not defined.
  3. A test that uses those variables and cannot say why it failed.
  4. Randomness with no seed — the one that produces no error at all.

The fourth is the dangerous one, because nothing tells you about it. The script runs, gives an answer, and gives a different answer to the next person.

NoteHint 2

The shape of the fixes:

  • Paths should be relative to the project, so "data/file.csv" works for anyone who has the repo. Absolute paths are a promise about someone else’s hard disk.
  • Uncomment the two definitions. They belong in the script because everything after them depends on them.
  • set.seed(2847) goes immediately before replicate(). Once, not inside the loop — a seed inside would make all thousand resamples identical.
TipSolution 
wind_data <- read.csv("data/uk_electricity.csv")

early <- wind_data$wind_twh[wind_data$year < 2015]
late <- wind_data$wind_twh[wind_data$year >= 2015]

result <- t.test(late, early)
result

set.seed(2847)
boot_means <- replicate(1000, mean(sample(late, replace = TRUE)))
quantile(boot_means, c(0.025, 0.975))

The bootstrap interval is 55.2 to 75.1 TWh, and now it is the same interval every time anyone runs it.

Three of the four faults announce themselves with an error. The fourth does not: unseeded randomness produces a working script whose numbers change quietly between runs. If you quote a bootstrap interval in your report and then rerun your analysis after a reviewer’s comment, the numbers will not match the ones you wrote, and you will not know why.

A fifth fault is not in this script and will be in yours: a library() call you never wrote down, because you loaded the package once, months ago, and your session has had it ever since. The only way to catch it is to restart R and run the script top to bottom, which is worth doing once a week and mandatory before you submit.

set.seed() is not about getting the “right” answer. It is about getting a checkable one — and everything in this module that separates evidence from assertion works the same way.

Exercise 1 continued: see the shape, not just the range

Run Exercise 1 first: this builds on its boot_means.

quantile(boot_means, c(0.025, 0.975)) threw away everything except two numbers. Plot boot_means itself and look at the shape those two numbers were cut from.

NoteHint 1

One column, mean_twh, on the x-axis. Which geom_*() shows the distribution of a single numeric column as bars?

NoteHint 2 
ggplot(boot_df, aes(x = mean_twh)) +
  geom_histogram(______)
TipSolution 
boot_df <- data.frame(mean_twh = boot_means)

ggplot(boot_df, aes(x = mean_twh)) +
  geom_histogram(bins = 30)

That histogram is what “55.2 to 75.1 TWh” is a summary of: the middle 95% of 1,000 resampled means. The two quantiles are two vertical slices through this shape, and the shape itself — where it’s tall, where it tapers — is information the two-number interval throws away.

One more thing to try: comment out set.seed(2847) in Exercise 1 above, rerun that block, then rerun this plot. The histogram shifts slightly and the interval changes with it — the shape is exactly as random as the two numbers were. Put the seed back before you move on; everything downstream in this page expects the seeded answer.

Extension: break your own script

Optional, and the highest-value hour of the week.

  1. Restart R — properly, so the workspace is empty — and run your own project script from the first line. Write down every error. That list is your reproducibility debt, and today is much cheaper than the night before the deadline.
  2. Put sessionInfo() at the end of your script and read what it prints. Which of those package versions would you need to record for someone to reproduce your numbers in three years?
  3. Find a number in your draft report and trace it back to the exact line of code that produced it. If you cannot, that number is not evidence yet. Do this for three numbers.
  4. Delete a variable your script depends on — rm(something) — and run again. Does the script fail loudly, or carry on and produce a plausible wrong answer? The second is the outcome to fear, and Exercise 3 is entirely about it.

Quick check

When you report a t-test result to a policy audience, what does the sentence need?

    1. The p-value
    1. The p-value and an effect size
    1. “The result was statistically significant”
    1. The effect size, a confidence interval, and a plain-language reading of what it means — with the p-value alongside

Set answer to "a", "b", "c" or "d".

Exercise 2: Say it in one sentence

You have the test from Exercise 1. Now produce the number a reader actually needs — the difference in mean annual wind generation between 2015–2025 and 2000–2014 — and then write the sentence.

NoteHint 1

The number is a subtraction between two means you can compute from elec — or read out of the t-test object you already have, whose estimate component holds both.

For the sentence: a reader who is numerate but not statistical needs to know which way, how much, how sure, and why it happened. In that order, and in one sentence.

NoteHint 2

The shape of it:

mean(late) - mean(early)

or, from the test object,

diff(rev(result$estimate))

For the confidence interval, result$conf.int gives 43.4 to 68.8. Round for the reader: nobody needs 43.40985.

TipSolution 
mean(late) - mean(early)

56.13 TWh. A model sentence:

UK wind generation averaged about 56 TWh a year more in 2015-2025 than in 2000-2014 (95% CI 43 to 69 TWh; Welch t = 9.42, p < 0.001) — roughly six times the earlier average, and equivalent to about a fifth of all UK electricity in 2025 — reflecting the large-scale deployment of onshore and offshore capacity over the period.

Take it apart:

  • Direction and size in real units. “56 TWh a year more”, not “significantly higher”.
  • Uncertainty, rounded. 43 to 69. The reader learns the estimate is good to within about ±23%, which is the honest precision.
  • A comparator. “Six times” and “a fifth of UK electricity” are what make 54 a big number rather than just a number. This is the refrain the whole module has been drilling: compared to what?
  • A mechanism. The last clause says why, which is what stops the sentence reading as a coincidence.

And one thing it does not do: claim that wind deployment caused the rise in a sense any test here established. It is an observed difference between two periods, and everything else about Britain changed too.

Common failure modes: a sentence that reports only p; a sentence that says “proves”; a sentence with the CI on a transformed scale a reader cannot picture.

Extension: the same result, three audiences

Optional.

  1. Rewrite your sentence for a journal reviewer — who wants the test, its assumptions, and the exact numbers — and then for a local councillor, who wants to know whether to approve a wind farm. Same evidence, and almost no words in common.
  2. Write the version a campaigner would publish. Keep every number true. Then mark which choices you made to get the effect — that list is what to look for when you read someone else’s.
  3. Find a results sentence in your own draft report and check it against the four ingredients. Most drafts are missing the comparator, because the author has looked at the number so long that it feels large on its own.

Quick check

AI tools are useful in research, and not uniformly. Which of these is AI most reliable at?

    1. Judging whether your sample size is adequate
    1. Finding a syntax error in your R code
    1. Deciding which statistical test fits your question
    1. Judging whether your effect size is large enough to matter

Set answer to "a", "b", "c" or "d".

Exercise 3: The confident wrong answer

You asked a language model: “Using the factors and elec data frames, how much CO₂ did UK biomass electricity emit in 2025 under the full supply chain scenario, in Mt?”

It replied with this code and the conclusion “0 Mt — biomass electricity is carbon neutral under this scenario.”

bio_factor <- sum(factors$co2_kg_per_mwh[
  factors$fuel == "Biomass" & factors$scenario == "with_supply_chain"
])

bio_2025 <- elec$bioenergy_twh[elec$year == 2025] * 1e6 * bio_factor / 1e9
bio_2025

The code runs. It produces no warning. The answer is wrong.

Find the fault, and print the correct figure.

NoteHint 1

Do not read the code looking for the mistake. Run it, then print every intermediate value and ask whether each is what you expected.

Start with bio_factor. It should be one number, around 1060. Print length(bio_factor) as well as its value — the length is the tell.

NoteHint 2

unique(factors$fuel) shows what is actually in that column.

When a subsetting condition matches nothing, R does not object: it returns a zero-length vector. Most functions then carry on. sum() of nothing is 0, length() of it is 0, and mean() of it is NaN — which at least looks wrong, whereas 0 looks like an answer.

This is why the fault survives: every step behaves exactly as documented.

TipSolution 
unique(factors$fuel)

bio_factor <- factors$co2_kg_per_mwh[
  factors$fuel == "biomass" & factors$scenario == "with_supply_chain"
]
bio_factor

elec$bioenergy_twh[elec$year == 2025] * 1e6 * bio_factor / 1e9

43.5 Mt, not zero.

The fault is a capital letter. factors$fuel contains "biomass"; the AI wrote "Biomass". The condition matched no rows, the subset was a zero-length vector, and sum() of a zero-length vector is 0 — not an error, but the documented identity element. Every step did exactly what it promised.

Three things make this the worst kind of bug:

  • It is silent. No error, no warning, nothing to grep for in a log.
  • The sum() hid it. Without it the multiplication would have produced numeric(0) and printed nothing, which at least looks wrong. Wrapping a subset in sum() converts “no data” into “zero”, and those are very different claims.
  • The answer was plausible. Zero is exactly what UK carbon accounting records for biomass. A reader who already believes the conclusion will not check the code, and the model’s fluent explanation makes the check feel unnecessary.

This is what to take into your report. The risk from a language model is not that it writes code that fails — you would notice. It is that it writes code that runs and quietly answers a slightly different question. The defence is the same as for your own work: print the intermediate values, check each against what you expected, and be most suspicious when the answer is one you wanted.

Extension: make no-match loud

Optional. ?stopifnot, ?match.arg, ?identical, ?all.equal.

  1. Rewrite the lookup so a no-match stops the script instead of returning zero. stopifnot(length(bio_factor) == 1) is one line and would have caught this. Add that check to your own project wherever you pull a single value out of a table.
  2. Find two other ways R turns “nothing” into a number that looks like an answer. sum() and prod() are the obvious pair; max() is the interesting one, because it warns and returns something absurd.
  3. Ask a language model the same question about this dataset and read its code before its answer. Did it check that its subset matched anything? Almost none do, and it is not because they cannot — it is because nobody in the training data bothered either.
  4. factors$fuel == "Biomass" is the case-sensitivity trap. Find the equivalent trap in your own data: a trailing space, an inconsistent capital, a name spelled two ways across two files.

Exercise 4: How plausible was it before you tested?

One last look at HolmesCo. Their gold assay is 95% accurate in both directions, and it returned 50 positives from 10,000 samples. Their press release calls these “50 confirmed gold-bearing sites”.

Suppose only 5 samples in 10,000 genuinely bear gold.

Print the positive predictive value: of the samples that test positive, what percentage really are gold-bearing?

NoteHint 1

Do not reach for Bayes’ theorem. Count people — or in this case, samples.

Imagine all 10,000 laid out. Five are gold-bearing; 9,995 are not. Send every one through the test and count the positives that come out of each pile. Then ask what fraction of the total positive pile came from the gold pile.

That reframing — natural frequencies rather than probabilities — makes this problem easy, and it is the reason the fallacy persists: the same arithmetic in probability notation defeats most people, including most doctors, most lawyers and most geologists.

NoteHint 2

The shape of it:

true_pos <- 5 * 0.95
false_pos <- (10000 - 5) * 0.05
100 * true_pos / (true_pos + false_pos)

Look at the two numbers before you divide. One is under 5, the other is about 500. The answer is decided by the sizes of the two piles, not by the accuracy of the test.

TipSolution 
true_pos <- 5 * 0.95
false_pos <- (10000 - 5) * 0.05

100 * true_pos / (true_pos + false_pos)

0.94%. About 4.75 real hits are hidden among roughly 500 false ones, so fewer than one positive in a hundred means anything. HolmesCo’s “50 confirmed sites” is a number their own test cannot support — and in fact 50 positives is fewer than the ~500 this test would throw up by chance alone, which is its own interesting question.

The general shape: when what you are looking for is rare, a good test still produces mostly false alarms, because the 5% error rate is applied to a vastly larger pile than the 95% success rate is.

This is the p-value problem in different clothing. A “significant” result for an implausible hypothesis is most likely a false positive, for exactly this arithmetic — which is why “how plausible was this before we tested?” has been a refrain since Week 1. It is not scepticism for its own sake. It is the denominator.

And it is the judgement an AI tool is least likely to volunteer. Ask one to analyse HolmesCo’s assay and it will usually compute what you asked for, fluently and correctly, without pointing out that the base rate makes the conclusion worthless. Noticing that the question was the wrong question remains your job.

Extension: what would make the test worth running?

Optional.

  1. How accurate would the test have to be for half the positives to be real, at this base rate? Solve it, and then say whether such a test could exist.
  2. Hold the test at 95% and raise the base rate instead. At what prevalence do half the positives become real? This is why screening is done on targeted populations rather than everybody, in medicine and in mineral exploration alike.
  3. HolmesCo found 50 positives where chance alone predicts about 500. Give two explanations, one innocent and one not. What would you ask to see in order to tell them apart?
  4. Write the sentence HolmesCo’s press release should have contained. Then ask why no consultancy writes it.

Save your work

The exercise that matters most this week is not on this page: restart R and run your own project script from the first line. Fix whatever breaks, commit the fixed version, and push via GitHub Desktop.

Then find one number in your draft report and trace it to the line of code that produced it. Do that for three numbers and you will have caught something.