---
title: "Computing in Statistics — Take-Home Midterm"
subtitle: "SOLUTION KEY | Fall 2026 | 100 Points"
author: "Instructor"
date: today
format:
  html:
    toc: true
    toc-depth: 2
    toc-title: "Contents"
    number-sections: true
    theme: cosmo
    highlight-style: github
    self-contained: true
  pdf:
    toc: true
    number-sections: true
    geometry: margin=1in
    fontsize: 11pt
execute:
  echo: true
  warning: false
  message: false
---

```{r setup, include=FALSE}
set.seed(42)
knitr::opts_chunk$set(fig.width = 7, fig.height = 4.2)
```

## Setup {.unnumbered}

```{r load-data}
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))
```

## Point Summary {.unnumbered}

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

---

# Exploring the Data *(36 points)*

## 1.1 Visual exploration *(16 pts)*

### Figure A — one histogram per variable

```{r}
#| fig-height: 7
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)
```

```{r}
sk <- vapply(df[plot_vars], function(z) mean((z - mean(z))^3)/sd(z)^3, 0)
round(sort(sk), 2)
```

**(a) Shape.** Roughly **symmetric**:
`r paste(names(sk)[abs(sk) < 0.35], collapse = ", ")`. Clearly
**right-skewed**: `r paste(names(sk)[sk > 1], collapse = ", ")`. They are not
alike, and the differences matter for which summary statistic is honest.

**(b) Mean versus median.** For BMI the mean is
`r round(mean(df$BMI), 2)` against a median of `r round(median(df$BMI), 2)` —
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.

```{r}
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))
```

`Poverty` is the cleanest answer: it is **capped at 5**, so
`r sum(df$Poverty == 5)` respondents
(`r round(100*mean(df$Poverty == 5), 1)`%) 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
(`r round(100*mean(df$DaysMentHlthBad == 0))`% and
`r round(100*mean(df$DaysPhysHlthBad == 0))`%).

::: {.callout-note collapse="true"}
## Marking note — Figure A (8 pts)
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.
:::

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

```{r}
#| fig-height: 5
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")
```

```{r}
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)
```

**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 **`r names(which.max(before))`**, skewness
`r round(max(before), 2)` both before and after: its median sits near the
bottom of its own range with a long run of outliers above.

::: {.callout-note collapse="true"}
## Marking note — Figure B (8 pts)
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 A summary function, applied to every measure *(12 pts)*

```{r}
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)
```

**Interpretation.**
`pct_above_mean` is `r round(s[["pct_above_mean"]], 1)`% — **far** from 50%. The
mean (`r round(s[["mean"]], 2)`) exceeds the median
(`r round(s[["median"]], 2)`) because BMI is bounded below but has a long upper
tail, with a maximum of `r round(s[["max"]], 1)`. 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.

::: {.callout-note collapse="true"}
## Marking note — part one (6 pts)
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]]`.
:::

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


```{r}
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")
```

**Interpretation.**
No — naming "the measure with the highest mean" is not meaningful. The means run
from `r round(min(summary_mat["mean", ]), 2)`
(`r names(which.min(summary_mat["mean", ]))`) to
`r round(max(summary_mat["mean", ]), 2)`
(`r names(which.max(summary_mat["mean", ]))`), 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.

::: {.callout-note collapse="true"}
## Marking note — part two (6 pts)
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 The highest and lowest value of each measure *(8 pts)*

```{r}
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")
```

```{r}
c(min_diastolic = min(dia), n_zero_diastolic = sum(dia == 0),
  pct = round(100 * mean(dia == 0), 2), min_pulse = min(pul))
```

**The unusual value.**

**(a)** `BPDiaAve` — diastolic blood pressure. Its minimum is **`r min(dia)`**.

**(b)** **`r sum(dia == 0)` respondents** out of `r nrow(df)`
(`r round(100*mean(dia == 0), 2)`%) 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:

```{r}
c(all_rows = mean(dia), zeros_removed = mean(dia[dia > 0]),
  shift = mean(dia[dia > 0]) - mean(dia))
```

The mean moves from `r round(mean(dia), 3)` to
`r round(mean(dia[dia>0]), 3)` mmHg, a shift of
`r round(mean(dia[dia>0]) - mean(dia), 3)` mmHg. Small — but that is a fact to
**establish**, not to assume.

Note `Pulse` has a minimum of `r min(pul)`, which is low but physiologically
possible in a resting adult; it is not the answer.

::: {.callout-note collapse="true"}
## Marking note (8 pts)
3 for a correct extremes table, 5 for the data check.

Full credit on the check needs: identifying `BPDiaAve == 0`, **counting** the
`r sum(dia == 0)` 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.
:::

---

# Classification & Group Summaries *(27 points)*

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

```{r}
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))
```

```{r}
# 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)
```

```{r}
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"))
```

**Interpretation.**
Check 3 matters because the two rules are **not mutually exclusive**: `r both`
respondents have systolic pressure at or above 130 *and* sleep under six hours.
The rules are priority-ordered, so all `r both` must come out as
`"Elevated BP"`. Had the sleep rule been tested first, exactly those `r both`
respondents would be mislabelled `"Short Sleep"` — about
`r round(100*both/nrow(df), 1)`% of the sample, and
`r round(100*both/sum(bp >= 130), 1)`% of everyone who should be in the
Elevated BP group.

::: {.callout-note collapse="true"}
## Marking note (14 pts)
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 `r both` 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 Profile composition by age decade *(13 pts)*

### Part one — overall counts and percentages

```{r}
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")
```

::: {.callout-note collapse="true"}
## Marking note — part one (3 pts)
`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.
:::

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


```{r}
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")
```

```{r}
tab
round(100 * prop.table(tab, 1), 1)
```

```{r}
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)
```

```{r}
pe <- 100 * prop.table(tab, 1)[, "Elevated BP"]
round(pe, 1)
```

**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 **`r round(pe[1], 1)`%** in the youngest
decade to **`r round(pe[length(pe)], 1)`%** in the oldest — a rise of
`r round(pe[length(pe)] - pe[1], 1)` percentage points, or roughly
`r round(pe[length(pe)]/pe[1], 1)`-fold.

**(c)** The mode discards everything except which category is largest. A
category can grow from `r round(pe[1], 1)`% to `r round(pe[length(pe)], 1)`% —
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.

::: {.callout-note collapse="true"}
## Marking note — parts two and three (10 pts)
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).
:::

---


# Resampling & Inference *(37 points)*

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

```{r}
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)
```

```{r}
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))
```

```{r}
tt <- t.test(bp[!active], bp[active]); tt
```

**Interpretation.**

**(a)** Inactive respondents averaged `r round(mean(bp[!active]), 2)` mmHg
against `r round(mean(bp[active]), 2)` for active respondents, a difference of
**`r round(T_obs, 2)` mmHg**. Of `r B` permutations, `r r_ext` reached that
magnitude, so the formula $r/B$ returns **`r p_perm`**. Stated correctly:
**among the respondents in this sample**, those who were not physically active
had higher mean systolic pressure.

Because `r r_ext` 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 = `r format(1/B, scientific = FALSE)`$.

**(b)** Yes, the two agree. Welch's t-test gives
`r format.pval(tt$p.value, digits = 3)`; 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 `r sum(!active)` and `r sum(active)` 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.

::: {.callout-note collapse="true"}
## Marking note (11 pts)
6 test, 5 interpretation.

::: {.callout-warning appearance="simple"}
**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 Is that difference what it appears to be? *(15 pts)*

```{r}
# Step 1 — is age related to the GROUPING?
round(tapply(df$Age, df$PhysActive, mean), 2)
round(100 * tapply(df$PhysActive == "No", df$AgeDecade, mean), 1)
```

```{r}
# Step 2 — is age related to the OUTCOME?
round(tapply(bp, df$AgeDecade, mean), 1)
round(cor(df$Age, bp), 3)
```

```{r}
# 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")
```

```{r}
# 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)))
```

**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
  `r round(mean(df$Age[df$PhysActive == "No"]), 1)` years against
  `r round(mean(df$Age[df$PhysActive == "Yes"]), 1)` for active ones, and the
  share who are inactive climbs steadily from
  `r round(100*mean(df$PhysActive[df$AgeDecade == decades[1]] == "No"), 1)`% in
  the youngest decade to
  `r round(100*mean(df$PhysActive[df$AgeDecade == decades[length(decades)]] == "No"), 1)`%
  in the oldest.
- *Condition (ii):* mean systolic pressure rises monotonically across decades,
  from `r round(tapply(bp, df$AgeDecade, mean)[1], 1)` to
  `r round(tapply(bp, df$AgeDecade, mean)[length(decades)], 1)` mmHg, with
  `cor(Age, BPSysAve)` = `r round(cor(df$Age, bp), 3)`.

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
**`r round(T_obs, 2)` mmHg**. Holding age fixed and recombining, it falls to
**`r round(adjusted, 2)` mmHg** — only
**`r round(100*adjusted/T_obs)`%** of the original, meaning roughly
**`r round(100*(1 - adjusted/T_obs))`% of the apparent effect was age**, not
activity.

**(c) Sign reversals.** In **`r sum(sign(within$difference) != sign(T_obs))` of
the `r length(decades)` 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.

::: {.callout-note collapse="true"}
## Marking note (15 pts)
**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).

::: {.callout-warning appearance="simple"}
**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 `r sum(sign(within$difference) != sign(T_obs))`
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 Bootstrap confidence interval for a difference in means *(11 pts)*

```{r}
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]])
```

```{r}
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
[`r round(ci_boot[[1]], 2)`, `r round(ci_boot[[2]], 2)`] mmHg, which
**`r if (ci_boot[[1]] <= 0 && 0 <= ci_boot[[2]]) "contains" else "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 `r format(signif(p_perm,3), scientific=FALSE)`$). 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.

::: {.callout-warning appearance="simple"}
**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
`r round(ci_boot[[1]], 1)` and `r round(ci_boot[[2]], 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.
:::

::: {.callout-note collapse="true"}
## Marking note (11 pts)
**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.
:::

### Step 6 — the bootstrap SE against the formula SE

```{r}
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)
```

**(c)** The bootstrap standard error is `r round(se_boot, 4)` mmHg; Welch's
formula gives `r round(se_ttest, 4)`. They agree to
`r round(100*abs(se_boot/se_ttest - 1), 1)`%. 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
[`r round(quantile(boot_diffs, .025), 2)`,
`r round(quantile(boot_diffs, .975), 2)`] on its own would repeat exactly the
error 3.2 was written to expose.

::: {.callout-note collapse="true"}
## Marking note — step 6 (2 of the 11)
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 {.unnumbered}

```{r}
sessionInfo()
```