Computing in Statistics — Take-Home Midterm (Programming Focus)

Fall 2026 | 100 Points | 30% of Class Grade

Author

Your Name

Published

October 7, 2026

Due Date: Oct 14, 2026 — before class

Instructions

Submit this .qmd file and its rendered output to .html. Then open it in a browser and print it to a .pdf file. Upload this pdf file to blackboard. All code must be visible (echo: true).
Any approach is allowed. Base R, the apply family (sapply(), lapply(), vapply(), mapply()), replicate(), ifelse(), case_when(), subset(), and the tidyverse (dplyr, tidyr, ggplot2, purrr, …) are all welcome. Use whichever tools produce clear, correct, reproducible code.
The code skeletons in each chunk are starting points, not requirements — replace them if you prefer a different style. You are graded on correctness, clarity, and interpretation, not on which construct you chose.
Provide a written interpretation after every computation that produces a result.
Every figure must be readable on its own. Give each plot an informative title, axis labels with units, and a legend where one is needed, and follow it with at least one sentence saying what it shows. A figure with no caption and no takeaway earns no credit, however pretty it is.
Never let a figure silently drop data. If you set axis limits or transform a scale, check that no points disappeared — R will usually warn you, and you should read those warnings.
Every figure you submit must be the output of a code chunk in your own document. Some questions show you a worked example of what the figure should look like. Those examples are watermarked, deliberately incomplete, and are not answers — a figure that is pasted in as an image, rather than drawn by your visible code, earns zero for that part.
Set set.seed(42) at the top and before every randomised block so your results are reproducible.

Permitted Resources & Academic Integrity

This is an open-resource, individual take-home exam. You are trusted to work independently, and the rules below exist so that trust is easy to honour.

You may use, without citing

All course materials — lecture notes, slides, lab handouts, the textbook, recorded lectures, the course website.
R’s own documentation — ?function, help(), package vignettes, and the messages R prints at you.
Your own prior work in this course.

You may use, but must cite

Any website: Stack Overflow, Cross Validated, package documentation, blog posts, tutorials, GitHub issues, video, anything else.
Any book or article outside the assigned textbook.
Any code you adapt from any external source, even a single line, even if you modified it.
Any factual claim or interpretation you took from elsewhere — including anything you look up about health measurements.

Searching the internet is expected, not frowned upon. Looking up how quantile() handles ties is exactly what a working statistician does. No marks are ever deducted for using a well-cited source. The violation is failing to cite, not consulting.

How to cite

Do both:

A comment at the point of use:

# Adapted from: https://stackoverflow.com/questions/XXXXXXX (accessed 2026-10-02)
# I used the approach of ordering the factor by median before plotting.

A ## References section at the end, listing each source once with a URL, the date accessed, and one line on what you used it for.

If a source shaped your interpretation rather than your code, cite it in the prose: “Following [3], I treated the median rather than the mean as the better summary for a skewed variable.”

Generative AI is not permitted

You may not use ChatGPT, Claude, Gemini, Copilot, or any other generative AI tool at any point while working on this exam. Specifically, you may not:

ask an AI tool to write, debug, explain, translate, comment, or review any code;
ask an AI tool to write, rewrite, expand, summarise, or polish any part of your interpretations;
paste any part of this exam — a question, your code, your output, or an error message — into any AI tool, including “just to check”;
leave AI autocomplete running in your editor. GitHub Copilot, the RStudio AI assistant, Cursor, and similar features generate code as you type. Turn them off before you begin, and say so in your declaration below.

Why this rule, when web search is allowed. Looking up how a function behaves leaves you with a fact you still have to apply. A generative tool hands over the applying — the loop, the interpretation, the sentence about what the histogram shows — which is precisely and only what this exam measures. Every question here asks you to make a modelling choice and justify it. An answer produced by a tool tells me nothing about whether you can make that choice, and tells you nothing either.

AI summaries in search results. Search engines now place AI-generated answers above the real results whether or not you asked for one. Seeing such a summary is not a violation — you did not seek it out. Relying on one is. Treat it as an unverifiable source with no author: scroll past it to the actual documentation or forum answer, confirm the claim there, and cite that. If you cannot find a real source for something an AI summary told you, do not use it.

Required declaration. End your submission, immediately above your ## References section, with this line:

I did not use any generative AI tool in preparing this submission, I disabled AI autocomplete in my editor while working on it, and I did not discuss this exam with anyone.

A submission without this line is incomplete and will be returned to you before it is graded.

Working with other people

This is an individual exam. You may not discuss it with anyone except me and the TA (Alex Goebel).

That prohibition covers classmates, students who have taken this course previously, tutors, friends, family, and anyone online. Specifically, you may not:

discuss this exam’s questions, concepts, or approaches with another person, in any medium — in person, by text, email, group chat, or a forum;
share, copy, view, or dictate code or written interpretations;
work through any question alongside someone else, even if no files change hands;
post any part of this exam anywhere, or ask about it on a public site.

Everything you submit must be your own work, written by you, in your own words.

The one exception is me and the TA (Alex Goebel). Questions about what a question is asking are welcome — bring them to office hours or email me. If an answer turns out to matter for everyone, I will email it to the whole class, so that nobody gains an advantage by having asked. Questions of the form “is my answer right?” are not appropriate on either channel.

If you get stuck

The routes open to you are the course materials, R’s documentation, a clearly-cited web source, office hours, and email to me — and failing all of those, a comment in your chunk explaining what you tried and where it broke. That comment earns partial credit. An AI-written answer, or an answer obtained by discussing the exam with someone, earns zero and is a violation of academic integrity, handled under the University’s policy rather than by me.

Practical requirements

Your .qmd must render without errors on a clean machine. Restart R and render from scratch before submitting. A submission that does not render loses the marks for anything that fails to run.
Do not hard-code results obtained elsewhere. Every number in your write-up must be produced by the code in your document.
If a question defeats you, submit your attempt with a comment on where you got stuck and what you tried. Partial credit is real; a blank chunk earns nothing.

Packages

# Load whatever you need, e.g.:
# library(ggplot2)
# library(dplyr)
# library(tidyr)
# library(knitr)

Dataset

NHANES — the US National Health and Nutrition Examination Survey. NHANES is a large, long-running study run by the Centers for Disease Control and Prevention, in which participants complete an interview and a physical examination: their height and weight are measured, their blood pressure is taken, and blood samples are analysed. That matters — most of the numbers below were measured by a technician, not self-reported.

You are given nhanes_wellness.csv, an extract of 5,397 adults aged 20–79 with no missing values. It was built from the NHANES R package, a teaching edition of the survey.

Important

Every conclusion you write should be a statement about the people in this dataset due to the fact that the survey does not sample the non-institutionalized civilian resident population of the United States evenly. It uses a complex survey design that deliberately oversamples certain subpopulations, including racial minorities, so the makeup of the raw data does not match the makeup of the country.

Defensible: “Among respondents in this sample, those who were physically active had lower mean systolic pressure.”
Not Defensible: “Physically active Americans have lower blood pressure.”

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

wellness_vars <- c("BMI", "Pulse", "BPSysAve", "BPDiaAve", "TotChol",
                   "DirectChol", "SleepHrsNight", "DaysMentHlthBad",
                   "DaysPhysHlthBad", "Poverty")

# Education is ordered; R will not know that unless you tell it
df$Education <- factor(df$Education,
  levels = c("8th Grade", "9 - 11th Grade", "High School",
             "Some College", "College Grad"))
df$AgeDecade <- factor(trimws(df$AgeDecade))

dim(df)

[1] 5397   20

head(df[, c("ID", "Age", "Gender", "BMI", "BPSysAve", "SleepHrsNight")], 3)

     ID Age Gender   BMI BPSysAve SleepHrsNight
1 51624  34   male 32.22      113             4
2 51624  34   male 32.22      113             4
3 51624  34   male 32.22      113             4

The ten wellness measures

These ten are used throughout. No health knowledge is required — each is described below, with a typical value so you can sanity-check your output.

Column	What it is	Units	Typical (median)
`BMI`	Body mass index: weight divided by height squared	kg/m²	27.9
`Pulse`	Resting heart rate	beats per minute	72
`BPSysAve`	Systolic blood pressure — the higher number, pressure when the heart beats	mmHg	118
`BPDiaAve`	Diastolic blood pressure — the lower number, pressure between beats	mmHg	71
`TotChol`	Total cholesterol in the blood	mmol/L	5.02
`DirectChol`	HDL cholesterol, the “good” kind	mmol/L	1.29
`SleepHrsNight`	Self-reported hours of sleep on a typical night	hours	7
`DaysMentHlthBad`	Days in the past 30 when mental health was not good	days (0–30)	0
`DaysPhysHlthBad`	Days in the past 30 when physical health was not good	days (0–30)	0
`Poverty`	Household income as a multiple of the federal poverty line. 1 = exactly at the line; capped at 5	ratio	3.14

Blood pressure is written as systolic “over” diastolic, e.g. 118/71.

Supporting columns

Column	Meaning
`ID`	Respondent identifier
`Age`, `AgeDecade`	Age in years (20–79); decade band (`20-29` … `70+`)
`Gender`	`female` / `male` as recorded by the survey
`Weight`, `Height`	Measured weight (kg) and height (cm) — plotted in Figure A
`Education`	Highest level completed, five ordered levels
`PhysActive`	Whether the respondent does moderate or vigorous recreational activity (`Yes`/`No`)
`SleepTrouble`	Whether they have told a doctor they have trouble sleeping
`HealthGen`	Self-rated general health: `Excellent`, `Vgood`, `Good`, `Fair`, `Poor`

Two things about these measures that will matter

They are on completely different scales. Cholesterol near 5, blood pressure near 120, sleep near 7. You cannot put them on one axis, or compare their means to each other, without doing something about that first. Figure B in Part 1.1 asks you to.

Some have hard boundaries. DaysMentHlthBad and DaysPhysHlthBad cannot exceed 30, and most people report 0. Poverty is capped at 5, so everyone above five times the poverty line is recorded as exactly 5. These caps shape the distributions you will plot, and they are properties of how the survey was recorded, not of the people.

stopifnot(all(wellness_vars %in% names(df)), !anyNA(df[, wellness_vars]))
length(wellness_vars)

[1] 10

Point Summary

One point = one percent of the midterm grade.
Part	Topic	Points	Weight
1	Exploring the Data	36	36%
2	Classification & Group Summaries	27	27%
3	Resampling & Inference	37	37%
Total		100	100%

How each question is laid out. Every question below has the same three headings, so you always know what is being asked of you:

What to do — the numbered steps. Follow them in order.
What to report — the specific numbers or objects that must appear in your output.
What to write — the questions your interpretation must answer. Answer each one explicitly.

1 Exploring the Data (36 points)

1.1 Visual exploration (16 pts)

Before computing a single summary statistic, look at the data. Produce two figures. Each needs a title, axis labels with units, and at least one sentence saying what it shows.

Figure A — one histogram per variable

What to do

Build a grid of histograms, one panel per numeric variable: Age, Weight, Height, and all ten measures in wellness_vars — thirteen panels.
Title each panel with its variable name.
Use free axis scales. The variables have very different ranges, and forcing a common axis makes most panels unreadable.

Base R: par(mfrow = c(3, 5)) followed by a for loop over the variable names, then par(op) to restore.

# YOUR CODE HERE — Figure A

What the example above shows, and what it does not. It shows the layout you are aiming at: thirteen panels, one per variable, each on its own scale. It is greyscale, watermarked, and its axis labels and units have been stripped out — your figure must have them. It is not an answer, and pasting it in earns zero.

What to write

Which distributions are roughly symmetric, and which have a long tail to the right? Name at least two of each.
Pick one variable. Add vertical lines at its mean and its median. Say which is larger, and explain why in terms of the shape you see.
Which distributions are shaped by a recording limit rather than by the people being measured? Name at least one and state what the limit is.

[Write your answer here.]

Figure B — a boxplot of all ten measures together

The ten measures are in different units, so you cannot put them on one axis as they stand. Standardise them first.

What to do

Convert each of the ten measures in wellness_vars to a z-score, $(x - \bar{x})/s$, so each has mean 0 and standard deviation 1. scale() does this in one call.
Reshape to long format: one row per (measured variable name, z-score) pair.
Draw a boxplot — geom_boxplot() in ggplot2, or boxplot() in base R — with the ten measures on one axis and the z-score on the other. coord_flip() makes the measure names readable.
Order the measures by spread, not alphabetically.
Do not set axis limits. If you do, some points will be dropped and R will warn you. Read the warning.

# YOUR CODE HERE — Figure B (boxplots)

Read this before you copy anything. The example above is ordered alphabetically. Step 4 requires yours to be ordered by spread. The example is therefore a figure that would lose marks if you submitted it — it is there to show you the form of the plot, not its content.

What to write

Does standardising change the shape of a distribution, or only its scale? Point to one specific measure in your boxplot that demonstrates your answer.

[Write your answer here.]

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

What to do — part one: write the function

Write a function with exactly this signature:

wellness_summary <- function(df, var) { ... }

Inside the function, extract the column named var from df. Refer to it programmatically — df[[var]] — not by typing the column name.
Return a named numeric vector with exactly these seven elements, in this order:

Element	Definition
`n`	number of non-missing observations
`mean`	arithmetic mean
`sd`	standard deviation
`median`	median
`min`	minimum
`max`	maximum
`pct_above_mean`	percentage of respondents whose value strictly exceeds the mean

wellness_summary <- function(df, var) {
  # YOUR CODE HERE
}

wellness_summary(df, "BMI")

NULL

What to do — part two: apply it to all ten

Call wellness_summary() on each of the ten measures in wellness_vars.
Assemble the results into a numeric matrix with 7 rows and 10 columns — one row per summary statistic, one column per measure.
Name the rows and columns, and display the matrix as a formatted table (Hint: knitr::kable(round(summary_mat, 2), caption = “Summary statistics by measure”))

A for loop, sapply(), or purrr::map_dfc() are all acceptable.

# YOUR CODE HERE — the 7 x 10 matrix

What to report

The seven summary statistics for "BMI", and the full 7 × 10 matrix.

What to write

Is pct_above_mean for BMI close to 50%, or far from it? Explain what that tells you about the shape of the BMI distribution, referring explicitly to the mean and the median. Connect your answer to the histogram in Figure A.
Look at the mean row of your matrix. Would it be meaningful to say which measure has “the highest mean”? Explain what stops the ten means from being comparable to one another, and name one thing you would have to do first for such a comparison to make sense.

[Write your answer here.]

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

What to do

For each of the ten measures in wellness_vars, find the respondent with the highest value and the respondent with the lowest value.
Collect the results in a data frame called extremes with exactly these five columns: measure, max_id, max_value, min_id, min_value.
Display it as a formatted table.

which.max() and which.min() return the position of the extreme value. Use that position to look up both the value and the respondent’s ID.

extremes <- data.frame(measure = wellness_vars, max_id = NA_integer_,
                       max_value = NA_real_, min_id = NA_integer_,
                       min_value = NA_real_)

# YOUR CODE HERE — fill in `extremes`

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	NA	NA	NA	NA
Pulse	NA	NA	NA	NA
BPSysAve	NA	NA	NA	NA
BPDiaAve	NA	NA	NA	NA
TotChol	NA	NA	NA	NA
DirectChol	NA	NA	NA	NA
SleepHrsNight	NA	NA	NA	NA
DaysMentHlthBad	NA	NA	NA	NA
DaysPhysHlthBad	NA	NA	NA	NA
Poverty	NA	NA	NA	NA

What to write — read this carefully

Look down the min_value and max_value columns and ask, for each one, “is this a value a real person could have?”

At least one of these numbers is not possible for a living adult. Find it. Then answer:

Which measure, and what is the impossible value?
How many respondents have it? (Count them — do not guess.)
What do you think it actually represents? Give the most likely explanation.

This is a data-checking exercise. Every dataset you will ever analyse contains something like this.

[Write your answer here.]

2 Classification & Group Summaries (27 points)

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

What to do

Write a function with exactly this signature:

classify_wellness <- function(sbp, sleep) { ... }

It takes one respondent’s systolic blood pressure and hours of sleep, and returns one character string, following these priority-ordered rules:

Returns	Condition
`"Elevated BP"`	`sbp` ≥ 130
`"Short Sleep"`	`sleep` < 6 — checked only if the first rule did not apply
`"Neither"`	neither condition holds

Apply your function to every respondent in df to produce a character vector of length nrow(df). A for loop over rows, orsapply(), will work.
Add the result to df as a new column called profile.

These are descriptive labels for this exercise, not diagnoses. 130 mmHg is the stage-1 hypertension threshold and under six hours is a common short-sleep cut-off, but a single reading does not diagnose anyone.

Important: the rules are priority-ordered, and a substantial number of respondents satisfy both conditions. Those must come out as "Elevated BP". Your function must get this right for a single respondent before you apply it to all of them.

classify_wellness <- function(sbp, sleep) {
  # YOUR CODE HERE — returns ONE string for ONE respondent
}

# Test your function on individual values before applying it to everyone.
# Fill in what each call should return, then check that it does:
classify_wellness(145, 7)    # expect "Elevated BP"

NULL

classify_wellness(145, 4)    # expect "Elevated BP"  (both rules -> first wins)

NULL

classify_wellness(110, 4)    # expect "Short Sleep"

NULL

classify_wellness(110, 8)    # expect "Neither"

NULL

profile <- character(nrow(df))   # placeholder so the blank template renders

# YOUR CODE HERE — apply classify_wellness() to every respondent

df[["profile"]] <- profile
head(df[, c("ID", "BPSysAve", "SleepHrsNight", "profile")], 10)

      ID BPSysAve SleepHrsNight profile
1  51624      113             4        
2  51624      113             4        
3  51624      113             4        
4  51630      112             8        
5  51647      118             8        
6  51647      118             8        
7  51647      118             8        
8  51654      111             7        
9  51656      104             5        
10 51657      134             4

What to report

Run these three checks and show the output:

every respondent received a label, and none is NA;
the count of respondents satisfying both rules;
How many respondents satisfy sbp ≥ 130 and sleep < 6? What are their wellness profile values?

# YOUR CODE HERE — the three checks

What to write

Why does check 3 matter? Describe what would go wrong — and how many respondents would be affected — if you had tested the sleep rule first.

[Write your answer here.]

2.2 Profile composition by age decade (13 pts)

What to do — part one: the overall counts

Count the respondents in each profile, and compute each count as a percentage of the total. Display both in a small table.

profile_labels <- c("Elevated BP", "Short Sleep", "Neither")

# YOUR CODE HERE — counts and percentages

When you build a table from profile, make sure the rows line up with the labels. table() sorts alphabetically, which is not the order of profile_labels. Setting factor(profile, levels = profile_labels) fixes it.

What to do — part two: the dominant profile in each decade

For each age decade, find which profile is most common.
Build a data frame with columns AgeDecade and dominant_profile, and display it.
State in a comment how your method breaks ties.

Your table should look like this (this is the correct answer — use it to check your code):

AgeDecade	dominant_profile
20-29	Neither
30-39	Neither
40-49	Neither
50-59	Neither
60-69	Neither
70+	Elevated BP

decade_profile <- data.frame(AgeDecade = levels(df$AgeDecade),
                             dominant_profile = NA_character_)

# YOUR CODE HERE — fill in `dominant_profile`

knitr::kable(decade_profile, caption = "Dominant Wellness Profile by Age Decade")

Dominant Wellness Profile by Age Decade
AgeDecade	dominant_profile
20-29	NA
30-39	NA
40-49	NA
50-59	NA
60-69	NA
70+	NA

What to do — part three: the figure

That table is nearly constant, and it hides almost everything. Build a 100%-stacked bar chart showing the percentage of each profile within every age decade.

The figure should have:

age decade on the horizontal axis, one bar per decade;
bars that each total 100%, divided into three coloured segments;
percentage on the vertical axis;
the three profiles distinguished by the colours given below.

Steps:

Build a table of counts: table(df$AgeDecade, df$profile).
Convert to row-wise proportions with prop.table(tab, margin = 1), so each decade sums to 1.
Convert to a data frame for plotting with as.data.frame(). You will get three columns; rename them AgeDecade, Profile and Freq.
Plot with geom_col(), multiplying Freq by 100 to get percentages.

Use this colour scheme so everyone’s figure is comparable:

scale_fill_manual(values = c("Elevated BP" = "#4C72B0",
                             "Short Sleep" = "#DD8452",
                             "Neither"     = "#BFBFBF"))

A skeleton of the plotting call:

ggplot(comp, aes(x = AgeDecade, y = 100 * Freq, fill = Profile)) +
  geom_col(width = 0.8) +
  scale_fill_manual(values = c("Elevated BP" = "#4C72B0",
                               "Short Sleep" = "#DD8452",
                               "Neither"     = "#BFBFBF")) +
  labs(title = "...", x = "...", y = "...", fill = NULL)

# YOUR CODE HERE — build `comp`, then draw the stacked bar chart

What the example above shows, and what it does not. It shows the form: one bar per decade, each reaching 100%, three colour bands in the scheme given above. Its title, axis labels and legend title have been removed — your figure must have all three, and as it stands this one would lose marks. It is not an answer, and pasting it in earns zero.

What to report

The overall counts and percentages from part one, and for each decade the percentage in each of the three profiles.

What to write

Your table of dominant profiles is almost the same in every row. Your figure is not. Describe the pattern the figure reveals that the table conceals.
Give the percentage in "Elevated BP" for the youngest decade and for the oldest. How large is the change?
Explain in one or two sentences why reporting only “the most common profile” was a poor way to summarise this relationship.

[Write your answer here.]

3 Resampling & Inference (37 points)

This part is worth 37 points. The three questions run in sequence on the same comparison: you will test it, discover that the test answered a narrower question than it appeared to, and then estimate how large the effect actually is. Answer them in order.

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

Test whether mean systolic blood pressure differs between respondents who are physically active and those who are not.

\[H_0: \mu_{\text{inactive}} = \mu_{\text{active}} \qquad H_a: \mu_{\text{inactive}} \neq \mu_{\text{active}}\]

What to do

Compute the observed statistic $T_{\text{obs}} = \bar{X}_{\text{inactive}} - \bar{X}_{\text{active}}$.
Set B <- 2000. Repeat B times: randomly shuffle the PhysActive labels across respondents, leaving the blood-pressure values where they are, and recompute the difference in group means. Store the results in perm_diffs.
Compute the two-sided p-value as $(r )/(B)$, where $r$ is the number of permuted differences with $|{\text{difference}}| \ge |T_{\text{obs}}|$.
Draw a histogram of perm_diffs with vertical lines at $T_{\text{obs}}$ and at $-T_{\text{obs}}$.
Run t.test() on the same comparison.

Report the permuted p-value. With B = 2,000 the smallest value you can obtain is $1/(B) \approx 0.0005$. If no permutation is as extreme as what you observed, the correct statement is $p \le 0.0005$.

set.seed(42)

bp     <- df$BPSysAve
active <- df$PhysActive == "Yes"

# Step 1: observed statistic
# YOUR CODE HERE for T_obs

# Step 2: permutation
B <- 2000
perm_diffs <- numeric(B)   # placeholder
# YOUR CODE HERE — fill `perm_diffs`

# Step 3: p-value
# YOUR CODE HERE

# Step 4: histogram
# YOUR CODE HERE

# Step 5: parametric comparison
t.test(bp[!active], bp[active])


    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

What to report

$T_{\text{obs}}$ in mmHg; the count $r$; your permutation p-value; and the p-value from t.test().

What to write

State your conclusion. Word it as a statement about the respondents in this sample, not about Americans in general — see the dataset note above.
Do the permutation test and the t-test agree? Which of the two relies on stronger assumptions, and what exactly are those assumptions?

[Write your answer here.]

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

In 3.1 you found that inactive respondents had higher mean systolic pressure, with a p-value small enough to look conclusive. This question asks whether that comparison means what it appears to mean.

Physically active and inactive respondents are not otherwise similar people. If some third variable differs between the two groups and is related to blood pressure, then part of what you measured in 3.1 could be due to that variable rather than to physical activity. A variable that does both of those things is a confounder.

Here you will check whether age is confounding the comparison.

What to do

Is age related to the grouping? Compute mean Age for active and inactive respondents separately. Then compute, for each level of AgeDecade, the percentage of respondents who are inactive.
Is age related to the outcome? Compute mean BPSysAve within each AgeDecade, and compute cor(Age, BPSysAve).
Hold age fixed. For each age decade separately, recompute the difference in mean BPSysAve between inactive and active respondents. Build a table with one row per decade and columns: decade, number of respondents, mean for the inactive, mean for the active, and the difference.
Combine the within-decade differences into a single number by taking a weighted average, weighting each decade’s difference by the number of respondents in that decade.
```
w <- table(df$AgeDecade) / nrow(df)
adjusted <- sum(w * within_decade_differences)
```
Compare that adjusted figure with the unadjusted difference from 3.1, and compute what percentage of the original difference remains.

# Step 1: is age related to the grouping?
# YOUR CODE HERE

# Step 2: is age related to the outcome?
# YOUR CODE HERE

# Steps 3-4: within-decade differences, then the weighted average
# YOUR CODE HERE

# Step 5: adjusted vs unadjusted
# YOUR CODE HERE

What to report

The two mean ages; the percentage inactive in each decade; mean blood pressure in each decade; cor(Age, BPSysAve); your full within-decade table; the weighted average; the unadjusted difference from 3.1; and the percentage of the unadjusted difference that survives adjustment.

What to write

A confounder must satisfy two conditions. State them, then show — using your own numbers from steps 1 and 2 — that age satisfies both here.
How much of the unadjusted difference survives once age is held fixed? Give the two numbers and the percentage change.
Look down the difference column of your within-decade table. In how many decades does the difference change sign compared with the overall result? What does that tell you about the claim you might have made after 3.1 alone?
Write the conclusion you would now defend about physical activity and blood pressure among the respondents in this sample. Be precise about what the data do and do not support.

A note on what you have and have not done. Holding age fixed by comparing within decades deals with age. It does nothing about any other variable that might differ between active and inactive respondents. Nothing in this exam establishes that physical activity causes a change in blood pressure, and your answer to (d) should not claim otherwise.

[Write your answer here.]

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

In 3.1 you tested whether mean blood pressure differs by physical activity. A p-value answers that yes-or-no question but says nothing about size. Here you will estimate how large the unadjusted difference plausibly is.

What to do

Compute the observed difference in means, $T_{\text{obs}} = \bar{X}_{\text{inactive}} - \bar{X}_{\text{active}}$.
Set B <- 2000. Repeat B times: draw nrow(df) rows of df with replacement, and recompute the difference in mean BPSysAve between the two activity groups within that resample. Store the results in boot_diffs.
Report the bootstrap standard error — the standard deviation of boot_diffs.
Report a 95% percentile confidence interval — the 2.5th and 97.5th percentiles of boot_diffs.
Draw a histogram of boot_diffs with vertical lines at $T_{\text{obs}}$ and at zero.
Extract the standard error that t.test() reported in 3.1 — it is t.test(bp[!active], bp[active])$stderr — and compare it with your bootstrap standard error.

Note what is being resampled here: whole rows, not the two groups separately. Each bootstrap replicate is a fresh draw of respondents from the dataset, and the group sizes within a replicate come out however they come out.

set.seed(42)

T_obs <- mean(bp[!active]) - mean(bp[active])

B <- 2000
boot_diffs <- numeric(B)   # placeholder

# YOUR CODE HERE — fill `boot_diffs`

# SE, 95% percentile CI, histogram
# YOUR CODE HERE

# Step 6: the two standard errors side by side
# YOUR CODE HERE

What to report

$T_{\text{obs}}$ in mmHg; the bootstrap standard error; the 95% percentile confidence interval; and the t.test() standard error from step 6.

What to write

Does the interval contain zero? State your conclusion, worded as a statement about the respondents in this sample.
Compare your conclusion with the permutation p-value from 3.1. Do the two agree?
How close is your bootstrap standard error to the one t.test() computed from a formula? What does that agreement demonstrate about what the bootstrap is doing?
This interval is for the unadjusted difference. In light of 3.2, what should a reader be told alongside it?

[Write your answer here.]

4 Declaration and References (required)

A submission missing either of these two sections is incomplete

It will be returned to you before it is graded. Fill both in before you render for the last time.

Academic Integrity Declaration

Copy the statement below into the blank beneath it, typed out in full, and add your name and the date. Typing it rather than ticking a box is the point: you are asserting it, not acknowledging it.

I did not use any generative AI tool in preparing this submission, I disabled AI autocomplete in my editor while working on it, and I did not discuss this exam with anyone other than the instructor or the TA.

Your declaration:

[Type the statement here.]

Name: [Your name]
Date: [Date]

If you cannot sign this truthfully, do not sign it. Come and talk to me instead. A conversation about what happened is recoverable; a false declaration is not.

References

List every source you consulted: websites, books or articles beyond the assigned textbook, and any code you adapted. One entry each, with a URL where there is one, the date you accessed it, and one line saying what you used it for.

Course materials, R’s own help files, and your own earlier work in this course do not need to be listed.

Format — delete the examples and replace them with your own:

Stack Overflow, “How to order a factor by median before plotting.” https://stackoverflow.com/questions/XXXXXXX — accessed 2026-10-09. Used for the reorder() call in Figure B.
R documentation, ?quantile. Consulted to check how ties are handled when computing the bootstrap percentile interval in 3.3.

[Your references here.]

If you used no outside sources, write that. “No external sources were consulted” is a complete and perfectly good References section. An empty heading is not — it reads as an oversight rather than a statement.

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     

loaded via a namespace (and not attached):
 [1] htmlwidgets_1.6.4 compiler_4.4.0    fastmap_1.2.0     cli_3.6.5        
 [5] tools_4.4.0       htmltools_0.5.9   rstudioapi_0.16.0 yaml_2.3.12      
 [9] rmarkdown_2.27    knitr_1.47        jsonlite_2.0.0    xfun_0.44        
[13] digest_0.6.39     rlang_1.1.7       evaluate_1.0.5