Week 6: Check Your Assumptions

When a significant result isn’t what it seems

Introduction

In the content session you learned how to run a t-test and read its output. But every statistical test rests on assumptions, and when those are violated the test can hand you a confident answer that is wrong.

Today’s data come from a HolmesCo borehole survey. HolmesCo measured temperatures at depth in two geological formations — the Whin Sill (a dolerite intrusion) and the Stainmore Formation (sedimentary) — to assess geothermal potential in County Durham.

By the end you will have produced a significant result, and then destroyed it, without touching a single number.

Goals:

  • Look at distributions before testing them
  • Run a t-test and read what it actually claims
  • Check normality rather than assuming it
  • Build the reflex: look, test, check, and be willing to go back

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: Look before you test

borehole has 70 rows and two columns: formation and temperature_c. Before running anything, find out what the two groups look like.

Draw one histogram per formation, then finish the block by printing the standard deviation of each group.

NoteHint 1

facet_wrap() splits a plot into panels by the values of a column, so you do not have to filter the data twice. It takes a formula: ~ column_name.

For the second part, you want one summary per group. Doing it by hand — subsetting each formation and calling sd() — works and is worth writing once. Then find the function that does it in a single call.

NoteHint 2

The shape of it:

ggplot(data, aes(x = value_column)) +
  geom_histogram(bins = 12) +
  facet_wrap(~ group_column)

tapply(data$value_column, data$group_column, sd)

tapply(values, groups, f) splits the first argument by the second and applies f to each piece, returning a named vector — which is why the output tells you which number belongs to which formation.

NoteHint 3 
ggplot(borehole, aes(x = ______)) +
  geom_histogram(bins = 12) +
  facet_wrap(~ ______)

tapply(borehole$temperature_c, borehole$______, ______)
TipSolution 
ggplot(borehole, aes(x = temperature_c)) +
  geom_histogram(bins = 12) +
  facet_wrap(~ formation)

tapply(borehole$temperature_c, borehole$formation, sd)

Both distributions are right-skewed — a long tail of high temperatures — and Stainmore is much the worse, with a few very hot boreholes pulling its mean upward. That is typical of geoscience measurements like temperature at depth, permeability and grain size, which are often roughly log-normal.

The spreads are 9.1 and 20.3: the Stainmore group varies more than twice as much. That is worth knowing but is not the problem. R’s t.test() uses Welch’s version by default, which handles unequal variances perfectly well.

The skew is the problem, and no default will save you from it.

Extension: four views of two groups

Optional. ?geom_boxplot, ?geom_violin, ?geom_density, ?scale_x_log10, ?geom_jitter.

  1. Show the same comparison without facet_wrap() — both groups in one panel. At least three geoms will do it. Which makes the difference in spread easiest to see, and which makes the difference in centre easiest?
  2. A boxplot of 35 points hides 35 points. Overlay the raw observations on it. Now count how many Stainmore boreholes are above 60 °C, and ask whether a summary of this group is meaningful at all.
  3. Put the x axis on a log scale and redraw the histograms. The shape changes completely — and you have not altered a single value. Keep that plot; Exercise 4 explains it.
  4. Which single chart would you put in front of a client who wanted a yes-or-no answer about geothermal potential? Which would you put in an appendix? They should not be the same chart.

Exercise 2: The significant result

Now do what HolmesCo did: run a two-sample t-test asking whether mean temperature differs between the formations.

End the block by printing the test.

NoteHint 1

R’s formula interface lets you say “this outcome, split by that group” in one expression. The tilde ~ means “described by”.

You do not need to split the data yourself — passing data = borehole lets R find both columns by name.

NoteHint 2

The shape of it:

t.test(outcome_column ~ grouping_column, data = df)

Print the result of that, rather than assigning it and printing one element. The output has four things worth reading: the t statistic, the degrees of freedom, the p-value, and the confidence interval for the difference in means.

The degrees of freedom will not be a whole number. That is Welch’s correction for unequal variances, and it is R’s default.

NoteHint 3 
raw_test <- t.test(______ ~ ______, data = borehole)
raw_test
TipSolution 
raw_test <- t.test(temperature_c ~ formation, data = borehole)
raw_test

p = 0.030. Significant at the conventional threshold, and HolmesCo would stop here.

Look at the confidence interval before you do: roughly 0.9 to 15.9 °C. The data are consistent with the Stainmore Formation being 16 °C hotter on average, and equally consistent with it being one degree hotter. A result can be “significant” and still be almost uninformative about the size of the thing you care about, which is why a p-value on its own is never a finding.

And there is a worse problem. You saw the skew in Exercise 1. The t-test assumes something about these data that you have not checked.

Extension: the test has options

Optional. ?t.test is the whole exercise — it is a longer help page than you expect, and almost every argument in it changes the answer.

  1. Run it again with var.equal = TRUE. That is Student’s original test, which pools the variances. The p-value barely moves here — construct group sizes where it would move a lot. (Hint: unequal variances hurt most when the group sizes are also unequal.)
  2. What does alternative = "greater" do to the p-value, and why is choosing it after seeing the data indefensible?
  3. The output reports degrees of freedom of about 48, not 68. Find out in the help page where the missing 20 went.
  4. t.test(x, y) and t.test(y ~ g) are two ways to call the same function. Reproduce your result using the other form, and check the sign of the estimate. Which group did R subtract from which, and how would you know without being told?

Exercise 3: Check the assumption

The t-test assumes the values within each group are roughly normally distributed. You have seen that they are not. Now demonstrate it, rather than asserting it.

A QQ plot puts your data against what a normal distribution would have produced. Points on the diagonal mean normal; a curve away from the line at one end means skew.

Finish by printing a formal test of normality for the Stainmore group.

NoteHint 1

Two of these steps you have done before. Subsetting one group’s values is the x[condition] pattern from Week 1, with a character comparison instead of a numeric one.

qqnorm() draws the plot and qqline() adds the reference line to the plot that already exists — the same “start it, then add to it” pattern as plot() and lines().

The third step is a test whose null hypothesis is the data are normal. Think about what a small p-value means when that is the null.

NoteHint 2

The shape of it:

group <- df$value[df$group_column == "GroupName"]

qqnorm(group)
qqline(group)

shapiro.test(group)

Careful with the name: the formation is recorded as Stainmore, and a comparison against a slightly different spelling returns no rows and no error.

NoteHint 3 
stainmore <- borehole$temperature_c[borehole$formation == "______"]

qqnorm(stainmore, main = "Stainmore, raw")
qqline(______)

shapiro.test(______)
TipSolution 
stainmore <- borehole$temperature_c[borehole$formation == "Stainmore"]

qqnorm(stainmore, main = "Stainmore, raw")
qqline(stainmore)

shapiro.test(stainmore)

p = 0.032. The null hypothesis of Shapiro-Wilk is that the data are normal, so a small p-value is evidence against normality — the opposite reading from the t-test you just ran, and a reliable source of confusion. Say it out loud each time: “small p, reject the null, the null was normality.”

The QQ plot shows the same thing without a threshold: the points curve away above the line at the top end, which is what a long right tail looks like.

Two cautions before you start applying this everywhere. First, the Whin Sill group passes (p = 0.28), so “the data are skewed” was really “one of the two groups is skewed”. Second, Shapiro-Wilk gets more powerful with sample size, so with 35 points it can miss real skew, and with 3,500 it will reject normality over deviations too small to matter. The QQ plot does not have that problem, which is a good argument for looking at one every time.

Extension: what does non-normal cost you?

Optional.

  1. Run Shapiro-Wilk on the Whin Sill group as well, and put its QQ plot beside Stainmore’s. ?par and its mfrow argument will get them side by side. Does the pair of plots change your view of Exercise 2?
  2. Generate 35 values from a genuinely normal distribution with rnorm(35, mean = 35, sd = 20) and QQ-plot them. Do it five times. How wiggly does a normal sample look at this size? That is the baseline you should be comparing against, and almost nobody establishes it before judging their own plot.
  3. Shapiro-Wilk on 3,500 values will reject normality for deviations far too small to affect a t-test. Demonstrate that with rnorm() plus a tiny amount of skew, and then decide which of the two tools — the test or the plot — you would rely on in a report.

Exercise 4: Transform and re-test

Taking logs pulls in a long right tail. If the values are roughly log-normal — as temperature at depth often is — the logged values will be roughly normal, and the t-test becomes entitled to its assumption.

Redo the test on logged temperatures and see what survives.

NoteHint 1

log() applies to a whole vector at once, so one line creates the new column. Assign it into the data frame with $ so that the formula interface can find it by name.

Then repeat Exercise 3’s check and Exercise 2’s test, with the new column in place of the old one. The only thing that changes is which column you name.

NoteHint 2

The shape of it:

df$log_value <- log(df$value)

shapiro.test(df$log_value[df$group == "..."])

t.test(log_value ~ group, data = df)

log() in R is the natural logarithm. Base 10 would work just as well for this purpose — the p-value is identical, because taking logs in a different base multiplies every value by a constant.

NoteHint 3 
borehole$log_temp <- ______(borehole$temperature_c)

shapiro.test(borehole$log_temp[borehole$formation == "Stainmore"])

t.test(______ ~ formation, data = borehole)
TipSolution 
borehole$log_temp <- log(borehole$temperature_c)

shapiro.test(borehole$log_temp[borehole$formation == "Stainmore"])

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

Shapiro-Wilk on the logged Stainmore values gives p = 0.21 — normality no longer rejected. The t-test then gives p = 0.22, nowhere near significant.

This is what happened to HolmesCo with their permeability survey: raw values gave p = 0.04 and a press release, logged values gave p = 0.23 and nothing. The original finding was driven by a handful of extreme observations in the tail, not by a difference between formations.

One honest complication, which matters more than the headline. The logged test is not the raw test made more reliable — it is a different question. A t-test on logs compares geometric means, so you are now asking whether the typical Stainmore borehole is hotter, rather than whether the average is. For skewed data that is usually the better question. But you have to say which one you answered, and switching to it after seeing the first p-value is a decision a reader is entitled to know about.

The lesson: a p-value is worth exactly as much as the assumptions behind it. Check first, and write down what you checked — including the checks you ran that did not change anything.

Extension: three more ways to handle a skewed tail

Optional. ?wilcox.test, ?sqrt, ?boxplot.stats.

  1. A Wilcoxon rank-sum test makes no normality assumption at all. Run it on the raw temperatures. You will get p = 0.09 — a third answer, from the same data, on the same question. Which of the three would you report, and could you defend that choice if you had chosen it before running any of them?
  2. Try a square-root transformation instead of a log. Does it fix the Shapiro-Wilk result? Transformations are not interchangeable, and the right one depends on the shape you are trying to fix.
  3. The alternative nobody should take but everyone considers: drop the extreme values and rerun the raw t-test. Do it, note the p-value, and then write down the sentence you would have to include in the methods section to describe what you did honestly. That sentence is why this is not an acceptable route.

Quick check

Across Exercises 2 and 4, did the p-value go up or down when you log-transformed?

Set answer to "increase" or "decrease".

Save your work

Copy the code you are most proud of into your week6.R file. Commit and push via GitHub Desktop. Write a commit message that describes what you did — not just “week 6”.