Computing in Statistics — Take-Home Midterm

SOLUTION KEY | Fall 2026 | 100 Points

Author

Instructor

Published

October 26, 2026

Setup

library(ggplot2)
theme_set(theme_minimal(base_size = 11))

df <- read.csv("nhanes_wellness.csv", stringsAsFactors = FALSE)

wellness_vars <- c("BMI", "Pulse", "BPSysAve", "BPDiaAve", "TotChol",
                   "DirectChol", "SleepHrsNight", "DaysMentHlthBad",
                   "DaysPhysHlthBad", "Poverty")
df$Education <- factor(df$Education,
  levels = c("8th Grade", "9 - 11th Grade", "High School",
             "Some College", "College Grad"))
df$AgeDecade <- factor(trimws(df$AgeDecade))

bp <- df$BPSysAve; dia <- df$BPDiaAve; pul <- df$Pulse
wt <- df$Weight;   ht  <- df$Height;   slp <- df$SleepHrsNight
active <- df$PhysActive == "Yes"

c(respondents = nrow(df), columns = ncol(df))
respondents     columns 
       5397          20

Point Summary

Part Topic Points
1 Exploring the Data 36
2 Classification & Group Summaries 27
3 Resampling & Inference 37
Total 100

1 Exploring the Data (36 points)

1.1 1.1 Visual exploration (16 pts)

1.1.1 Figure A — one histogram per variable

plot_vars <- c("Age", "Weight", "Height", wellness_vars)
op <- par(mfrow = c(3, 5), mar = c(3.4, 3.2, 2.2, 1), mgp = c(2, .7, 0))
for (v in plot_vars)
  hist(df[[v]], breaks = 25, col = "#4C72B0", border = "white",
       main = v, xlab = NULL, ylab = "Count", cex.main = .95)
par(op)

sk <- vapply(df[plot_vars], function(z) mean((z - mean(z))^3)/sd(z)^3, 0)
round(sort(sk), 2)
       BPDiaAve   SleepHrsNight         Poverty          Height             Age 
          -0.62           -0.30           -0.15            0.06            0.15 
          Pulse         TotChol        BPSysAve          Weight      DirectChol 
           0.48            0.68            0.97            1.02            1.05 
            BMI DaysMentHlthBad DaysPhysHlthBad 
           1.23            2.32            2.72

(a) Shape. Roughly symmetric: Age, Height, SleepHrsNight, Poverty. Clearly right-skewed: Weight, BMI, DirectChol, DaysMentHlthBad, DaysPhysHlthBad. They are not alike, and the differences matter for which summary statistic is honest.

(b) Mean versus median. For BMI the mean is 28.97 against a median of 27.9 — the mean is larger, dragged up by the long right tail visible in its panel.

(c) Shaped by a recording limit. Three qualify; students need name only one.

c(DaysMentHlthBad_zeros = sum(df$DaysMentHlthBad == 0),
  DaysPhysHlthBad_zeros = sum(df$DaysPhysHlthBad == 0),
  Poverty_at_cap        = sum(df$Poverty == 5),
  pct_Poverty_at_cap    = round(100 * mean(df$Poverty == 5), 1))
DaysMentHlthBad_zeros DaysPhysHlthBad_zeros        Poverty_at_cap 
               3098.0                3579.0                1424.0 
   pct_Poverty_at_cap 
                 26.4

Poverty is the cleanest answer: it is capped at 5, so 1424 respondents (26.4%) sit at exactly that value. That spike is not a cluster of people with identical incomes — it is everyone above five times the poverty line recorded as the same number. The two “bad days” measures are bounded at 30 with large spikes at zero (57% and 66%).

The grid must cover all thirteen variables with free scales. Part (c) is the discriminator: the student must name a variable whose shape comes from how it was recorded, not from the people measured. Accept Poverty’s cap or either “bad days” ceiling. A student who describes the Poverty spike as “many people with the same income” has the observation but not the insight — half credit.

1.1.2 Figure B — a boxplot of all ten measures, standardised

Z <- scale(df[, wellness_vars])
long <- data.frame(measure = rep(wellness_vars, each = nrow(df)),
                   z = as.vector(Z))
ord <- names(sort(tapply(long$z, long$measure,
                         function(q) diff(quantile(q, c(.25, .75))))))
long$measure <- factor(long$measure, levels = ord)

ggplot(long, aes(measure, z)) +
  geom_boxplot(fill = "grey90", outlier.size = .4, outlier.alpha = .25) +
  coord_flip() +
  labs(title = "All ten wellness measures, standardised to z-scores",
       subtitle = "Ordered by interquartile range. No points clipped.",
       x = NULL, y = "Standard deviations from the mean")

before <- vapply(df[wellness_vars], function(z) mean((z-mean(z))^3)/sd(z)^3, 0)
after  <- apply(Z, 2, function(z) mean((z-mean(z))^3)/sd(z)^3)
round(rbind(before, after), 3)
         BMI Pulse BPSysAve BPDiaAve TotChol DirectChol SleepHrsNight
before 1.232 0.483    0.973   -0.625   0.684      1.051        -0.301
after  1.232 0.483    0.973   -0.625   0.684      1.051        -0.301
       DaysMentHlthBad DaysPhysHlthBad Poverty
before           2.317           2.719   -0.15
after            2.317           2.719   -0.15

Answer. Standardising changes the scale, not the shape. Every measure now has mean 0 and SD 1 and they share an axis — but subtracting a constant and dividing by a constant cannot shorten a long tail, and the skewness row above is identical before and after to three decimals. The clearest example in the boxplot is DaysPhysHlthBad, skewness 2.72 both before and after: its median sits near the bottom of its own range with a long run of outliers above.

The question said draw a boxplot. A density plot or histogram grid does not answer it. Standardising is required — plotting the ten raw measures together is the common error and produces a meaningless figure. Reward ordering by spread.

Watch for silent data loss. A student who adds ylim() to tidy the axis drops BMI’s tail and gets a ggplot warning saying so. The paper tells them to read warnings. Deduct if points vanished with no comment; give credit to anyone who notices and either keeps them or states the clipping.

1.2 1.2 A summary function, applied to every measure (12 pts)

wellness_summary <- function(df, var) {
  x <- df[[var]]; x <- x[!is.na(x)]
  m <- mean(x)
  c(n = length(x), mean = m, sd = sd(x), median = median(x),
    min = min(x), max = max(x),
    pct_above_mean = 100 * length(x[x > m]) / length(x))
}
s <- wellness_summary(df, "BMI")
round(s, 3)
             n           mean             sd         median            min 
      5397.000         28.968          6.679         27.900         15.020 
           max pct_above_mean 
        81.250         43.617

Interpretation. pct_above_mean is 43.6% — far from 50%. The mean (28.97) exceeds the median (27.9) because BMI is bounded below but has a long upper tail, with a maximum of 81.2. A minority of high values pulls the mean away from the bulk of the data, so fewer than half the sample sits above it. Under a symmetric distribution the figure would be near 50%; the gap from 50% is the skew, and it is the same skew visible in BMI’s panel in Figure A.

5 for a correct function, 2 for the interpretation. Reserve the interpretation marks for linking pct_above_mean to the mean–median gap and back to Figure A. Bugs: >= instead of strict >; returning a list rather than a named numeric vector; hard-coding the column name instead of using df[[var]].

1.2.1 Part two — the 7 × 10 matrix across all ten measures

summary_mat <- sapply(wellness_vars, function(v) wellness_summary(df, v))
stopifnot(is.matrix(summary_mat), dim(summary_mat) == c(7, 10))
knitr::kable(round(summary_mat, 2), caption = "Summary statistics by measure")
Summary statistics by measure
BMI Pulse BPSysAve BPDiaAve TotChol DirectChol SleepHrsNight DaysMentHlthBad DaysPhysHlthBad Poverty
n 5397.00 5397.00 5397.00 5397.00 5397.00 5397.00 5397.00 5397.00 5397.00 5397.00
mean 28.97 72.49 120.10 70.78 5.09 1.37 6.87 4.08 3.29 3.04
sd 6.68 11.67 16.00 11.93 1.06 0.42 1.32 7.79 7.39 1.65
median 27.90 72.00 118.00 71.00 5.02 1.29 7.00 0.00 0.00 3.14
min 15.02 40.00 78.00 0.00 1.53 0.39 2.00 0.00 0.00 0.00
max 81.25 128.00 226.00 116.00 13.65 4.03 12.00 30.00 30.00 5.00
pct_above_mean 43.62 45.41 44.12 52.84 46.45 41.56 63.63 23.92 20.86 51.36

Interpretation. No — naming “the measure with the highest mean” is not meaningful. The means run from 1.37 (DirectChol) to 120.1 (BPSysAve), but in different units: mmHg, mmol/L, hours, days, and a dimensionless ratio. Asking whether 120 millimetres of mercury exceeds 5 millimoles per litre has no answer.

To compare them you must remove the units — standardise, as Figure B does — or compare each measure only against itself, as Part 5 does.

1 for the matrix, 2 for the units point. A student who answers “BPSysAve has the highest mean” without qualification has missed the question entirely.

1.3 1.3 The highest and lowest value of each measure (8 pts)

extremes <- data.frame(measure = wellness_vars, max_id = NA_integer_,
                       max_value = NA_real_, min_id = NA_integer_,
                       min_value = NA_real_)
for (i in seq_along(wellness_vars)) {
  x <- df[[wellness_vars[i]]]; h <- which.max(x); l <- which.min(x)
  extremes$max_id[i] <- df$ID[h]; extremes$max_value[i] <- x[h]
  extremes$min_id[i] <- df$ID[l]; extremes$min_value[i] <- x[l]
}
knitr::kable(extremes, caption = "Highest and Lowest Value of Each Measure")
Highest and Lowest Value of Each Measure
measure max_id max_value min_id min_value
BMI 52315 81.25 61853 15.02
Pulse 66589 128.00 54092 40.00
BPSysAve 55311 226.00 58821 78.00
BPDiaAve 59487 116.00 53529 0.00
TotChol 62146 13.65 67117 1.53
DirectChol 67487 4.03 60332 0.39
SleepHrsNight 51952 12.00 53243 2.00
DaysMentHlthBad 51702 30.00 51654 0.00
DaysPhysHlthBad 51723 30.00 51624 0.00
Poverty 51647 5.00 51982 0.00 
c(min_diastolic = min(dia), n_zero_diastolic = sum(dia == 0),
  pct = round(100 * mean(dia == 0), 2), min_pulse = min(pul))
   min_diastolic n_zero_diastolic              pct        min_pulse 
            0.00            13.00             0.24            40.00

The unusual value.

(a) BPDiaAve — diastolic blood pressure. Its minimum is 0.

(b) 13 respondents out of 5397 (0.24%) have a recorded diastolic pressure of zero.

(c) A living adult cannot have zero diastolic pressure — it would mean arterial pressure falling to nothing between heartbeats. These are measurement or coding artifacts: most likely a failed automated reading, or a missing value encoded as 0 rather than NA.

(d) Not deletion-on-sight. The sequence is: count them (done above); check the source documentation for whether 0 is a missing-value code for this instrument; recode to NA; re-run any analysis that uses the variable; and report that you did so, with the count. Deleting rows quietly is the one unacceptable option, because it changes results with no record.

(e) Yes — BPDiaAve is plotted in Figure A, standardised in Figure B, and supplies one column of the 1.2 matrix. The effect on that column:

c(all_rows = mean(dia), zeros_removed = mean(dia[dia > 0]),
  shift = mean(dia[dia > 0]) - mean(dia))
     all_rows zeros_removed         shift 
   70.7763572    70.9472511     0.1708939

The mean moves from 70.776 to 70.947 mmHg, a shift of 0.171 mmHg. Small — but that is a fact to establish, not to assume.

Note Pulse has a minimum of 40, which is low but physiologically possible in a resting adult; it is not the answer.

3 for a correct extremes table, 5 for the data check.

Full credit on the check needs: identifying BPDiaAve == 0, counting the 13 rather than saying “a few”, calling them artifacts rather than real observations, and proposing to investigate and recode rather than drop.

Part (e) is a bonus discriminator. A student who traces the value back into their own Figure A, Figure B or 1.2 matrix and quantifies the shift has done genuine analysis hygiene. Do not require it; reward it generously where present.

A student who flags Pulse = r min(pul) instead has been appropriately suspicious but picked the wrong target — partial credit, and note in feedback that a resting pulse that low is unusual but real.

2 Classification & Group Summaries (27 points)

2.1 2.1 A function that assigns a wellness profile (14 pts)

classify_wellness <- function(sbp, sleep) {
  if (sbp >= 130) {
    "Elevated BP"          # first rule wins even when sleep < 6 as well
  } else if (sleep < 6) {
    "Short Sleep"
  } else {
    "Neither"
  }
}

c(classify_wellness(145, 7), classify_wellness(145, 4),
  classify_wellness(110, 4), classify_wellness(110, 8))
[1] "Elevated BP" "Elevated BP" "Short Sleep" "Neither"    
# mapply() walks the two vectors in parallel, calling the scalar function once
# per respondent. A for loop over seq_len(nrow(df)) is equally correct.
profile <- mapply(classify_wellness, bp, slp)
df$profile <- profile
head(df[, c("ID", "BPSysAve", "SleepHrsNight", "profile")], 10)
      ID BPSysAve SleepHrsNight     profile
1  51624      113             4 Short Sleep
2  51624      113             4 Short Sleep
3  51624      113             4 Short Sleep
4  51630      112             8     Neither
5  51647      118             8     Neither
6  51647      118             8     Neither
7  51647      118             8     Neither
8  51654      111             7     Neither
9  51656      104             5 Short Sleep
10 51657      134             4 Elevated BP
both <- sum(bp >= 130 & slp < 6)
c(no_NA = !anyNA(df$profile),
  all_three_labels_only = all(df$profile %in%
      c("Elevated BP", "Short Sleep", "Neither")),
  satisfy_both_rules = both,
  all_both_labelled_ElevatedBP = all(df$profile[bp >= 130 & slp < 6] ==
                                     "Elevated BP"))
                       no_NA        all_three_labels_only 
                           1                            1 
          satisfy_both_rules all_both_labelled_ElevatedBP 
                         146                            1

Interpretation. Check 3 matters because the two rules are not mutually exclusive: 146 respondents have systolic pressure at or above 130 and sleep under six hours. The rules are priority-ordered, so all 146 must come out as "Elevated BP". Had the sleep rule been tested first, exactly those 146 respondents would be mislabelled "Short Sleep" — about 2.7% of the sample, and 11.9% of everyone who should be in the Elevated BP group.

5 for a correct scalar function, 1 for applying it to every respondent, 2 for the checks and the explanation.

The function must take one respondent’s values and return one string. The paper is explicit about this. A student who writes a vectorised ifelse() over the whole column has not answered as asked — deduct the function marks, though the classification logic may be perfectly correct and should still earn the logic marks.

The 146 overlapping respondents make a wrong rule order detectable. Run check 4 against their code rather than eyeballing the first ten printed rows, which look identical either way.

2.2 2.2 Profile composition by age decade (13 pts)

2.2.1 Part one — overall counts and percentages

profile_labels <- c("Elevated BP", "Short Sleep", "Neither")
pf <- factor(df$profile, levels = profile_labels)
results <- data.frame(
  profile        = profile_labels,
  count          = as.integer(table(pf)),
  proportion_pct = 100 * as.numeric(table(pf)) / nrow(df))
stopifnot(sum(results$count) == nrow(df))
knitr::kable(results, digits = 2, caption = "Wellness Profile Summary")
Wellness Profile Summary
profile count proportion_pct
Elevated BP 1226 22.72
Short Sleep 517 9.58
Neither 3654 67.70

factor(levels = profile_labels) matters: plain table() sorts alphabetically (Elevated BP, Neither, Short Sleep) and silently misaligns the counts against hard-coded labels. Check the rows match the numbers, not just that the numbers look plausible.

2.2.2 Parts two and three — dominant profile, then the stacked bar

tab <- table(df$AgeDecade, pf)
decade_profile <- data.frame(
  AgeDecade = rownames(tab),
  dominant_profile = colnames(tab)[apply(tab, 1, which.max)])
knitr::kable(decade_profile, caption = "Dominant Wellness Profile by Age Decade")
Dominant Wellness Profile by Age Decade
AgeDecade dominant_profile
20-29 Neither
30-39 Neither
40-49 Neither
50-59 Neither
60-69 Neither
70+ Elevated BP 
tab
       pf
        Elevated BP Short Sleep Neither
  20-29          81         114     829
  30-39         112          97     809
  40-49         208         133     774
  50-59         295          94     658
  60-69         308          47     398
  70+           222          32     186
round(100 * prop.table(tab, 1), 1)
       pf
        Elevated BP Short Sleep Neither
  20-29         7.9        11.1    81.0
  30-39        11.0         9.5    79.5
  40-49        18.7        11.9    69.4
  50-59        28.2         9.0    62.8
  60-69        40.9         6.2    52.9
  70+          50.5         7.3    42.3
comp <- as.data.frame(prop.table(tab, 1))
names(comp) <- c("AgeDecade", "Profile", "Freq")
ggplot(comp, aes(AgeDecade, 100 * Freq, fill = Profile)) +
  geom_col(width = .8) +
  scale_fill_manual(values = c("Elevated BP" = "#4C72B0",
                               "Short Sleep" = "#DD8452",
                               "Neither"     = "#BFBFBF")) +
  labs(title = "Wellness profile composition by age decade",
       x = "Age decade", y = "Share of respondents (%)", fill = NULL)

pe <- 100 * prop.table(tab, 1)[, "Elevated BP"]
round(pe, 1)
20-29 30-39 40-49 50-59 60-69   70+ 
  7.9  11.0  18.7  28.2  40.9  50.5

Interpretation.

(a) The table reads "Neither" in five of six rows and reveals nothing. The figure shows the "Elevated BP" band widening steadily and monotonically across every decade while "Neither" shrinks to match — a clean, strong, entirely orderly pattern that the dominant-profile column conceals completely.

(b) "Elevated BP" runs from 7.9% in the youngest decade to 50.5% in the oldest — a rise of 42.5 percentage points, or roughly 6.4-fold.

(c) The mode discards everything except which category is largest. A category can grow from 7.9% to 50.5% — becoming the single clearest pattern in the dataset — and still lose five rows out of six. Reporting only the winner threw away the entire finding.

1 table, 3 figure, 2 interpretation.

The table was given in the paper, so there is no credit for reproducing it — the 2 points are for producing it with correct code and stating the tie rule (which.max() returns the first maximum, so ties break by column order).

The figure must be 100%-stacked with bars summing to 100, using the supplied colours. Marks here are for the figure and for articulating what the mode threw away.

This is the clearest “a summary statistic can hide the whole pattern” example in the exam. A student who produces the figure but writes only “Elevated BP increases with age” has described it without answering (c).

3 Resampling & Inference (37 points)

3.1 3.1 Permutation test: blood pressure by physical activity (11 pts)

set.seed(42)
T_obs <- mean(bp[!active]) - mean(bp[active])
B <- 2000; perm_diffs <- numeric(B)
for (b in 1:B) { g <- sample(active)
  perm_diffs[b] <- mean(bp[!g]) - mean(bp[g]) }
r_ext  <- sum(abs(perm_diffs) >= abs(T_obs))
p_perm <- r_ext / B                      # the paper's formula, r/B
c(mean_inactive = mean(bp[!active]), mean_active = mean(bp[active]),
  T_obs = T_obs, r = r_ext, p_permutation = p_perm, resolution = 1/B)
mean_inactive   mean_active         T_obs             r p_permutation 
   121.365874    119.090000      2.275874      0.000000      0.000000 
   resolution 
     0.000500 
hist(perm_diffs, breaks = 40, col = "grey85", border = "white",
     main = "Permutation null: difference in mean systolic BP",
     xlab = "Inactive minus active (mmHg)")
abline(v = c(T_obs, -T_obs), col = "firebrick", lwd = 2, lty = c(1, 2))

tt <- t.test(bp[!active], bp[active]); tt

    Welch Two Sample t-test

data:  bp[!active] and bp[active]
t = 5.1417, df = 4852.7, p-value = 2.829e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.408123 3.143625
sample estimates:
mean of x mean of y 
 121.3659  119.0900

Interpretation.

(a) Inactive respondents averaged 121.37 mmHg against 119.09 for active respondents, a difference of 2.28 mmHg. Of 2000 permutations, 0 reached that magnitude, so the formula \(r/B\) returns 0. Stated correctly: among the respondents in this sample, those who were not physically active had higher mean systolic pressure.

Because 0 permutations were as extreme, \(r/B\) evaluates to exactly zero. Zero is not a defensible p-value — no finite simulation can establish that a result is impossible. The honest report is \(p < 1/B = 0.0005\).

(b) Yes, the two agree. Welch’s t-test gives 2.83e-07; both reject decisively.

The t-test assumes more. The permutation test needs only that the activity labels are exchangeable under \(H_0\). The t-test additionally requires the sampling distribution of the difference in means to be approximately normal. With 2397 and 3000 respondents per group the CLT makes that extra assumption nearly free here — which is why the two agree. At ten per group they could easily diverge.

6 test, 5 interpretation.

Instructor: a conflict in the paper. Step 3 asks for \(r/B\), which here is exactly 0, while the callout immediately below it says the smallest obtainable value is \(1/(B+1) \approx 0.0005\) and that \(p = 0\) must never be reported. Students will notice. Accept all three of: 0 with a sentence explaining it means \(p < 1/B\); \(p < 0.0005\); and \((r+1)/(B+1)\). Award full credit to anyone who flags the tension — that is the statistically literate response. Consider reconciling the two before the paper goes out.

The trap that remains: the conclusion must be about this sample. “Physically active Americans have lower blood pressure” over-claims; the paper warns twice and offers a worked example of correct wording — deduct.

The interpretation marks are mostly for the assumptions contrast in (b), not for observing that two p-values are both small.

Note: whether activity causes lower pressure is taken up directly in 3.2, which asks students to check age as a confounder. Here they are asked only whether the groups differ in this sample — a student who jumps ahead and raises confounding at 3.1 should be credited, not penalised.

3.2 3.2 Is that difference what it appears to be? (15 pts)

# Step 1 — is age related to the GROUPING?
round(tapply(df$Age, df$PhysActive, mean), 2)
   No   Yes 
48.43 43.79 
round(100 * tapply(df$PhysActive == "No", df$AgeDecade, mean), 1)
20-29 30-39 40-49 50-59 60-69   70+ 
 33.8  39.7  45.2  47.1  49.1  63.6 
# Step 2 — is age related to the OUTCOME?
round(tapply(bp, df$AgeDecade, mean), 1)
20-29 30-39 40-49 50-59 60-69   70+ 
113.2 115.5 118.6 123.4 126.7 131.7 
round(cor(df$Age, bp), 3)
[1] 0.364
# Steps 3-4 — hold age fixed, then recombine
decades <- levels(df$AgeDecade)
within <- data.frame(decade = decades, n = NA_integer_,
                     inactive = NA_real_, activ = NA_real_, difference = NA_real_)
for (i in seq_along(decades)) {
  s <- df[df$AgeDecade == decades[i], ]; a <- s$PhysActive == "Yes"
  within$n[i]          <- nrow(s)
  within$inactive[i]   <- mean(s$BPSysAve[!a])
  within$activ[i]      <- mean(s$BPSysAve[a])
  within$difference[i] <- mean(s$BPSysAve[!a]) - mean(s$BPSysAve[a])
}
knitr::kable(within, digits = 2,
             caption = "Inactive minus active systolic BP, within each age decade")
Inactive minus active systolic BP, within each age decade
decade n inactive activ difference
20-29 1024 111.76 113.88 -2.12
30-39 1018 114.78 115.89 -1.11
40-49 1115 118.96 118.36 0.60
50-59 1047 126.45 120.61 5.84
60-69 753 127.36 126.08 1.28
70+ 440 130.18 134.22 -4.05 
# Step 5 — adjusted vs unadjusted
w        <- within$n / sum(within$n)
adjusted <- sum(w * within$difference)
c(unadjusted = T_obs, adjusted = adjusted,
  pct_remaining = 100 * adjusted / T_obs,
  pct_removed = 100 * (1 - adjusted / T_obs),
  n_decades_reversing = sum(sign(within$difference) != sign(T_obs)))
         unadjusted            adjusted       pct_remaining         pct_removed 
          2.2758740           0.4955471          21.7739240          78.2260760 
n_decades_reversing 
          3.0000000

Interpretation.

(a) The two conditions. A confounder must be (i) associated with the grouping variable, and (ii) associated with the outcome — while not being a consequence of the grouping. Age satisfies both:

  • Condition (i): inactive respondents average 48.4 years against 43.8 for active ones, and the share who are inactive climbs steadily from 33.8% in the youngest decade to 63.6% in the oldest.
  • Condition (ii): mean systolic pressure rises monotonically across decades, from 113.2 to 131.7 mmHg, with cor(Age, BPSysAve) = 0.364.

Older people are both less likely to be active and have higher blood pressure. So some of what 3.1 attributed to activity is really age.

(b) How much survives. The unadjusted difference was 2.28 mmHg. Holding age fixed and recombining, it falls to 0.5 mmHg — only 22% of the original, meaning roughly 78% of the apparent effect was age, not activity.

(c) Sign reversals. In 3 of the 6 decades the difference points the opposite way from the overall result — in those decades the active respondents had the higher mean pressure. A single pooled number had concealed that the within-group pattern is not even consistent in direction. Any claim made after 3.1 alone would have been describing a pattern that does not hold in half the age bands it supposedly covers.

(d) A defensible conclusion. Among the respondents in this sample, inactive respondents had higher mean systolic pressure than active ones, but most of that difference is accounted for by age: inactive respondents are older on average, and blood pressure rises with age. After holding age fixed, the remaining difference is small and inconsistent in direction across age decades. Nothing here establishes that physical activity changes anyone’s blood pressure, and adjusting for age does nothing about any other variable that differs between the two groups.

7 for the two association checks and the within-decade table; 4 for the weighted recombination and the adjusted-vs-unadjusted comparison; 4 for the written parts (a)–(d).

Instructor: a part was dropped. The revised paper no longer asks why a small p-value offers no protection against confounding — the old part (d). That was the deepest question in the exam. If you want it back, re-add it to the paper and restore the answer below. What follows is kept for your own use in class.

Why the p-value gave no protection. A p-value answers exactly one question: could chance alone produce a difference this large? It is silent on a difference in what. Confounding is a property of the comparison, not of the sampling, so no p-value can detect it, and more data shrinks the p-value while leaving the confounding untouched. Do not accept “the p-value was wrong” — it was correct for the question it answered. Do not accept “p-values can be misleading” without naming which question goes unanswered.

(a) must name both conditions and attach the student’s own numbers to each. Naming them from memory without evidence from steps 1 and 2 is half credit.

(c) rewards actually reading the table. The 3 sign reversals are visible only if the student looks down the difference column; many will report the weighted average and stop.

On (d): the conclusion must retreat from 3.1’s claim without over-correcting into “there is no relationship.” The adjusted difference is small, not demonstrably zero. Deduct for either over-claim.

Accept any sensible weighting in step 4 — equal weights per decade, or weights by decade size as here. Students who weight by decade size and say why have done slightly better. Accept also a student who stratifies by a coarser age split (say under/over 50) provided they justify it; the shrinkage will be smaller but the qualitative finding survives.

Do not require students to notice that adjusting for age leaves other variables unaddressed — but reward it warmly where present, since it is the habit that makes the whole lesson stick.

3.3 3.3 Bootstrap confidence interval for a difference in means (11 pts)

set.seed(42)
B <- 2000; boot_diffs <- numeric(B)
for (b in 1:B) {
  i <- sample.int(nrow(df), nrow(df), replace = TRUE)
  boot_diffs[b] <- mean(bp[i][!active[i]]) - mean(bp[i][active[i]])
}
se_boot <- sd(boot_diffs); ci_boot <- quantile(boot_diffs, c(.025, .975))
c(T_obs = T_obs, bootstrap_SE = se_boot,
  lower = ci_boot[[1]], upper = ci_boot[[2]])
       T_obs bootstrap_SE        lower        upper 
   2.2758740    0.4439318    1.3884186    3.1301089 
hist(boot_diffs, breaks = 40, col = "grey85", border = "white",
     main = "Bootstrap distribution: difference in mean systolic BP",
     xlab = "Inactive minus active (mmHg)")
abline(v = T_obs, col = "firebrick", lwd = 2)
abline(v = 0, col = "grey40", lty = 3)

Interpretation.

(a) The 95% percentile interval is [1.39, 3.13] mmHg, which excludes zero. Among the respondents in this sample, inactive respondents had higher mean systolic pressure.

(b) This agrees with the permutation p-value from 3.1 (\(p \le 0\)). It should: a 95% interval that excludes zero and a two-sided test that rejects at the 5% level are close to the same statement made two ways — both ask whether zero is a plausible value for the difference, one by interval and one by tail probability.

Instructor: two parts were dropped from the paper. The revised 3.3 no longer asks what the interval tells a reader that the p-value cannot, nor what the bootstrap distribution represents. Both answers are kept here for class use; neither may be graded.

What the interval adds. A magnitude: the difference is plausibly between 1.4 and 3.1 mmHg, a scale a reader can judge against clinical relevance. The p-value says only “not zero”, and would have said that just as emphatically for 0.2 mmHg at this sample size.

What varies between resamples. Which respondents happen to be drawn. Each replicate is a fresh sample of the same size, drawn with replacement — the same kind of variation that would arise from collecting a new sample. The spread of boot_diffs therefore estimates the sampling variability of the difference, which is exactly what a standard error is.

5 for a correct bootstrap with SE and percentile interval, 2 for parts (a) and (b), 2 for the SE comparison in (c), 2 for the caveat the paper labels (e).

The implementation must resample whole rows, not the two groups separately to fixed sizes. Percentile, basic and normal-approximation intervals are all acceptable if the student names which one they used.

The discriminator is now the final part. A student who reports the interval without connecting it back to 3.2 has not understood what 3.2 established. The answer wanted is that resampling stands in for repeated sampling, so the spread of the replicates estimates sampling variability. Part (c) is the second discriminator: many students will say the interval “confirms” the p-value without noticing it adds a magnitude the p-value never had.

3.3.1 Step 6 — the bootstrap SE against the formula SE

se_boot  <- sd(boot_diffs)
se_ttest <- t.test(bp[!active], bp[active])$stderr
c(bootstrap = se_boot, t_test = se_ttest, ratio = se_boot / se_ttest)
bootstrap    t_test     ratio 
0.4439318 0.4426279 1.0029459

(c) The bootstrap standard error is 0.4439 mmHg; Welch’s formula gives 0.4426. They agree to 0.3%. That is the point of the bootstrap: it reproduces by brute force the spread the formula derives algebraically — and it keeps working for statistics that have no such formula.

(e) — the paper labels this part (e) although it follows (c); there is no (d) in the revised question. This interval is for the unadjusted difference. After 3.2 we know roughly 78% of it is age, so a reader must be told that the age-adjusted difference is far smaller and not distinguishable from zero. Reporting [1.39, 3.13] on its own would repeat exactly the error 3.2 was written to expose.

The comparison must be made, not merely printed. A student who reports both numbers without noticing they agree, or who does not connect (e) back to the confounding result, loses these two points. Part (e) is the one that separates a student who understood 3.2 from one who executed it.

Session Info

sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS 26.6.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_4.0.2

loaded via a namespace (and not attached):
 [1] vctrs_0.7.2        cli_3.6.5          knitr_1.47         rlang_1.1.7       
 [5] xfun_0.44          generics_0.1.3     S7_0.2.1           jsonlite_2.0.0    
 [9] labeling_0.4.3     glue_1.8.0         htmltools_0.5.9    scales_1.4.0      
[13] rmarkdown_2.27     grid_4.4.0         tibble_3.2.1       evaluate_1.0.5    
[17] fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5    compiler_4.4.0    
[21] dplyr_1.2.0        RColorBrewer_1.1-3 pkgconfig_2.0.3    htmlwidgets_1.6.4 
[25] rstudioapi_0.16.0  farver_2.1.2       digest_0.6.39      R6_2.6.1          
[29] tidyselect_1.2.1   dichromat_2.0-0.1  pillar_1.11.1      magrittr_2.0.3    
[33] withr_3.0.2        tools_4.4.0        gtable_0.3.6