Week 6: First Tests

Applying hypothesis tests to real data

Introduction

The content session took the borehole data apart one assumption at a time. This session is about putting the whole thing back together — running the workflow as a single piece of work, reading what it actually says, and then doing it on your own project data.

Today’s goals:

Run a complete t-test workflow in one block, unprompted
Read a confidence interval, including on a transformed scale
Write a results sentence a policy reader could use
See why a significant p-value is the start of the argument
Build a null world by shuffling labels, and see that it agrees with the t-test
Apply all of it to your group’s data

How these exercises work

The code boxes start empty; the numbered comments say what each step should produce. Hint 1 restates the goal, Hint 2 sketches the shape, Hint 3 gives you the code with the key pieces blanked, and the solution gives you the lot. Try before you open them.

Two things worth knowing:

The boxes share one R session. Anything you create in Exercise 1 is still there in Exercise 2.
The last thing your block prints is what gets marked, so variable names are yours to choose.

From this week on the blocks get longer. You have done every one of these steps individually — the new skill is assembling them without being walked through, because that is what analysing your own data will feel like.

Every exercise has an extension underneath it.

Stuck?

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: The whole workflow, one block

You did these four steps separately in the content session. Do them together now, as you would in a real script.

The question: do borehole temperatures differ between the Whin Sill and the Stainmore Formation?

Hint 2

The shape of it:

ggplot(df, aes(x = group, y = value)) + geom_violin() + geom_point()

shapiro.test(df$value[df$group == "A"])
shapiro.test(df$value[df$group == "B"])

df$log_value <- log(df$value)
t.test(log_value ~ group, data = df)

Or replace the last two lines with wilcox.test(value ~ group, data = df) and skip the transformation entirely.

Hint 3

ggplot(borehole, aes(x = ______, y = ______)) +
  geom_violin() +
  geom_point()

shapiro.test(borehole$temperature_c[borehole$formation == "Whin Sill"])
shapiro.test(borehole$temperature_c[borehole$formation == "______"])

# Stainmore fails the normality check, so:
borehole$log_temp <- ______(borehole$temperature_c)

t.test(______ ~ formation, data = borehole)

Solution

ggplot(borehole, aes(x = formation, y = temperature_c)) +
  geom_violin() +
  geom_point()

shapiro.test(borehole$temperature_c[borehole$formation == "Whin Sill"])
shapiro.test(borehole$temperature_c[borehole$formation == "Stainmore"])

# Stainmore is right-skewed (Shapiro-Wilk p = 0.03), and temperature at
# depth is plausibly log-normal, so test on the log scale.
borehole$log_temp <- log(borehole$temperature_c)

t.test(log_temp ~ formation, data = borehole)

The plot shows Stainmore sitting higher and spreading much further, with a long upper tail of a few very hot boreholes. Shapiro-Wilk confirms what that implies: Whin Sill passes at p = 0.28, Stainmore fails at p = 0.03. Logging fixes it, and the t-test then gives p = 0.22.

Look at what the comment in Step 3 is doing. It is not explaining the code — anyone can see that log() takes a logarithm. It records the decision, and it is the only part of this block that a reader could disagree with. Comments that narrate code are noise; comments that justify choices are the audit trail.

Extension: the order matters

Optional.

Run the workflow in the wrong order: test first, then check, then transform only if you did not like the answer. You will get the same three numbers. Write one sentence explaining why the sequence changes what those numbers are worth, even though the arithmetic is identical.
Wrap the whole workflow in a function of two arguments — a data frame and a column name — so you could point it at your project data tomorrow. ?function, and note that referring to a column whose name is in a variable needs df[[name]] rather than df$name.
Have your function print a warning when a group fails Shapiro-Wilk, rather than leaving you to notice. ?warning. Then ask whether you would actually read that warning, or scroll past it.

Exercise 2: What the interval says

The p-value told you the difference is not distinguishable from zero. It did not tell you how large a difference the data would still be consistent with — and that is usually the question a client is asking.

Your test reported a confidence interval on the log scale, where the numbers are meaningless to a reader. Undo the transformation so the interval describes a ratio of temperatures.

Finish by printing the back-transformed interval: two numbers.

Hint 2

The shape of it:

test_object$conf.int

exp(test_object$conf.int)

exp() applies to both ends at once, which is what you want: the transformation is monotonic, so the bounds stay bounds and stay in order.

Solution

log_test <- t.test(log_temp ~ formation, data = borehole)

log_test$conf.int

exp(log_test$conf.int)

On the log scale the interval is −0.09 to 0.40, which tells a reader nothing. Exponentiated it becomes 0.91 to 1.49: the Stainmore Formation’s typical temperature is between 9% cooler and 49% warmer than the Whin Sill’s.

Now the result is usable. It contains 1, so the data cannot establish a difference — but the upper end is 1.49, so they are also perfectly consistent with Stainmore being half again as hot. A client deciding whether to drill needs that second sentence far more than the p-value.

This is why “no significant difference” is a dangerous phrase. It is often heard as “the two are the same”. What it means here is “35 boreholes per group were not enough to tell”.

Extension: intervals are the interesting output

Optional. ?t.test for conf.level, ?names, ?str.

Recompute the interval at 99% confidence, then at 80%. Which is wider, and why is a wider interval sometimes the more useful one to report?
names(log_test) lists everything the object holds. Find the degrees of freedom and the standard error, and confirm by hand that the interval is the estimate plus or minus about two standard errors.
Work out how many boreholes per group you would need for the interval to exclude a 20% difference, assuming the spread stays as it is. ?power.t.test does this in one call. Then ask what it would have cost HolmesCo to collect them — this is exactly the calculation that should happen before a survey, and almost never does.

Exercise 3: Write the sentence

A result that nobody can read is not a result. Write one sentence reporting this analysis for a policy audience — someone numerate but not statistical, who has to make a decision.

It needs four things: what was compared, the size and direction of the difference, the uncertainty, and what follows. There is no code to write; type your sentence as a comment.

Extension: three bad sentences

Optional.

Each of these is defensible line by line and misleading overall. For each, name the trick in a few words, then rewrite it honestly.

“There was no significant difference in borehole temperature between the two formations (p = 0.22), confirming that the formations are geothermally equivalent.”
“Stainmore boreholes were 8.4 °C hotter on average than Whin Sill boreholes, a difference of over 30%.”
“Initial analysis showed a significant difference between formations (p = 0.03).”

The third one is the hardest, because nothing in it is false. Work out what a reader would have to be told for it not to mislead, and note that the fix is a sentence the author would rather not write.

Exercise 4: Significance is not size

With enough data, a difference too small to matter becomes “significant”. This is the single most consequential thing to understand about p-values, and it is easiest to see by building it.

Simulate two groups whose true means differ by 0.5 °C against a spread of 5 °C — a difference nobody would act on — and watch the p-value as the sample grows.

Hint 2

The shape of it:

set.seed(7042)

for (n in c(10, 30, 100, 500, 2000)) {
  a <- rnorm(n, mean = 25, sd = 5)
  b <- rnorm(n, mean = 25.5, sd = 5)
  p <- t.test(a, b)$p.value
  cat("n =", n, " p =", round(p, 4), "\n")
}

cat() rather than print() inside a loop, because print() labels every line with an index you do not want. And note that a loop’s last expression is not returned — you need a separate final line for the value being checked.

Hint 3

set.seed(______)

p_last <- NA
for (n in c(10, 30, 100, 500, ______)) {
  a <- rnorm(n, mean = 25, sd = 5)
  b <- rnorm(n, mean = ______, sd = 5)
  p_last <- t.test(a, b)$______
  cat("n =", n, " p =", round(p_last, 4), "\n")
}

p_last

Solution

set.seed(7042)

p_last <- NA
for (n in c(10, 30, 100, 500, 2000)) {
  a <- rnorm(n, mean = 25, sd = 5)
  b <- rnorm(n, mean = 25.5, sd = 5)
  p_last <- t.test(a, b)$p.value
  cat("n =", n, " p =", round(p_last, 4), "\n")
}

p_last

n = 10     p = 0.9625
n = 30     p = 0.1469
n = 100    p = 0.4313
n = 500    p = 0.0179
n = 2000   p = 0.0034

At n = 10 the difference is undetectable. At n = 2000 it is “highly significant” — and it is the same 0.5 °C difference throughout, worth nothing to anybody. Significance measures how well you can see an effect, not how much the effect matters.

Do not miss the untidy part. n = 100 gives a larger p-value than n = 30. That is not a mistake in the code; it is what sampling noise looks like. A p-value is a property of the sample you happened to draw, not a measurement of the world, and a single study’s p-value could easily have come out on the other side of 0.05 by luck alone.

Next week gives you the tool that fixes this: an effect size, which says how big the difference is and does not grow more impressive as you collect more data.

Extension: how unstable is a p-value?

Optional. ?replicate, ?power.t.test.

Run the n = 30 comparison 1,000 times with different seeds and histogram the p-values. What fraction fall below 0.05? Compare that with power.t.test(n = 30, delta = 0.5, sd = 5).
Now do the same with a true difference of zero. The histogram should be flat — every p-value equally likely. Convince yourself of that, because it is the cleanest statement of what a p-value is, and it explains exactly why running five tests and reporting the best one is cheating.
Add a Cohen’s d — the difference in means divided by the pooled standard deviation — to your original loop. Watch it stay put while the p-value collapses. That is the whole of next week in one column of numbers.

Exercise 5: Build the null world

In the content session the formation labels were shuffled a thousand times to build a world where formation cannot matter, and the p-value turned out to be nothing more than the share of that world at least as extreme as the real data. That was on raw temperatures. Exercise 1 settled on the log scale — so do it again there, and see whether the null world agrees with the t-test you trusted.

The five steps: measure δ, build a null world, drop δ into it, count, decide.

Hint 2

The shape of it:

log_temp <- log(df$value)
delta <- mean(log_temp[df$group == "A"]) - mean(log_temp[df$group == "B"])

null_deltas <- replicate(1000, {
  shuffled <- sample(df$group)
  mean(log_temp[shuffled == "A"]) - mean(log_temp[shuffled == "B"])
})

mean(abs(null_deltas) >= abs(delta))

mean() of a TRUE/FALSE vector is the fraction that are TRUE.

Hint 3

set.seed(7042)

log_temp <- ______(borehole$temperature_c)
delta <- mean(log_temp[borehole$formation == "Stainmore"]) -
  mean(log_temp[borehole$formation == "Whin Sill"])

null_deltas <- replicate(______, {
  shuffled <- ______(borehole$formation)
  mean(log_temp[shuffled == "Stainmore"]) -
    mean(log_temp[shuffled == "Whin Sill"])
})

mean(______(null_deltas) >= ______(delta))

Solution

set.seed(7042)

log_temp <- log(borehole$temperature_c)
delta <- mean(log_temp[borehole$formation == "Stainmore"]) -
  mean(log_temp[borehole$formation == "Whin Sill"])

null_deltas <- replicate(1000, {
  shuffled <- sample(borehole$formation)
  mean(log_temp[shuffled == "Stainmore"]) -
    mean(log_temp[shuffled == "Whin Sill"])
})

mean(abs(null_deltas) >= abs(delta))

δ is 0.152 on the log scale, and 220 of the 1,000 shuffled worlds produce a difference at least that large: p = 0.22. The logged t-test in Exercise 1 said 0.217. No formula, no normality assumption, no degrees of freedom — and the same answer.

That agreement is the whole lesson. Every test you will meet asks the same question: how surprising would this be if nothing were going on? The t-test answers it with a curve that someone worked out on paper; you answered it by building the world where nothing is going on and looking. When the two disagree, the curve’s assumptions are the first suspect.

Notice what the null world does not tell you: how likely it is that the formations really differ. That is a different question — the base-rate one — and it needs something the data alone cannot give you: how plausible the difference was before you tested.

Extension: δ can be anything

Optional. ?median, ?wilcox.test.

Swap the difference in means for the difference in medians, on raw temperatures. No textbook test is built for exactly that, and you do not need one — the null world works for any δ you can compute. How does your p-value compare with wilcox.test(temperature_c ~ formation, data = borehole)?
Run your Exercise 5 code again without the seed, five times. How much does the p-value wobble with 1,000 shuffles? With 10,000? How many would you need before the wobble stopped mattering for a decision at 0.05?
Draw the null world: a histogram of null_deltas with delta marked by geom_vline(), both signs. Shade or colour the tails if you can. That picture is the p-value, and it is a better figure for a report than the number on its own.

Now: your project data

The rest of the session is your group’s. Run the Exercise 1 workflow on your own data, and keep the order: look, check, decide, test, interpret.

Load and visualize. Points, violins or histograms by group. Anything surprising in the picture is more important than anything in the test.
Check the assumptions. Normal within each group? Skewed? How many observations per group — and are the groups the same size?
Decide, and write down why. Transform, or use a rank-based test, or proceed and note the caveat. All three can be right; only silence is wrong.
Test, and read the interval, not just the p-value.
Write it up. One paragraph, using Exercise 3’s four ingredients.

If your project has more than two groups, resist the urge to test every pair — that is next week’s topic and there is a trap in it. For today, pick the two groups your question is actually about, and say why those two.

If your data are not ready yet

You are not stuck. Do one of these instead, and all three are genuinely useful to the project:

Write the analysis plan as comments: which columns, which groups, which test, and what result would change your mind. Running it later becomes a half-hour job.
Generate fake data with the shape you expect using rnorm(), and run your whole workflow on it. If the code breaks, better now. If the result looks implausible even when you invented the data, your plan needs rethinking.
Take the borehole data and ask a question we did not: is the variance different between formations? ?var.test.

Save your work

Copy your analysis code into your week6.R file. Commit and push via GitHub Desktop. Include both the walkthrough and whatever project analysis you got to — and if you wrote a comment explaining a choice, make sure it survives the copy.