--- title: "Computing in Statistics — Take-Home Midterm (Programming Focus)" subtitle: "Fall 2026 | 100 Points | 30% of Class Grade" author: "Your Name" date: "Oct 7, 2026" format: html: toc: true toc-depth: 2 toc-title: "Contents" number-sections: true theme: cosmo highlight-style: github code-fold: false 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) ``` ::: {.content-visible when-format="html"} ::: {style="font-size:1.9em; font-weight:800; text-align:center; color:#b3001b; margin:1.1em 0 0.9em; letter-spacing:-0.01em;"} Due Date: Oct 14, 2026 — before class ::: ::: ::: {.content-visible when-format="pdf"} \begin{center} {\LARGE\bfseries Due Date: Oct 14, 2026 --- before class} \end{center} ::: --- ## Instructions {.unnumbered} - 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 {.unnumbered} 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 {.unnumbered} - **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** {.unnumbered} - **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 {.unnumbered} Do both: 1. **A comment at the point of use:** ```r # Adapted from: https://stackoverflow.com/questions/XXXXXXX (accessed 2026-10-02) # I used the approach of ordering the factor by median before plotting. ``` 2. **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** {.unnumbered} **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 {.unnumbered} ::: {.callout-important} ## 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 {.unnumbered} 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 {.unnumbered} - 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 {.unnumbered} ```{r packages} # Load whatever you need, e.g.: # library(ggplot2) # library(dplyr) # library(tidyr) # library(knitr) ``` ## Dataset {.unnumbered} **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. ::: {.callout-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."* ::: ```{r load-data} 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) head(df[, c("ID", "Age", "Gender", "BMI", "BPSysAve", "SleepHrsNight")], 3) ``` ### The ten wellness measures {.unnumbered} 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 {.unnumbered} | 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` | ::: {.callout-note} ## 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. ::: ```{r check-data} stopifnot(all(wellness_vars %in% names(df)), !anyNA(df[, wellness_vars])) length(wellness_vars) ``` ## Point Summary {.unnumbered} | Part | Topic | Points | Weight | |---|---|:---:|:---:| | 1 | Exploring the Data | 36 | 36% | | 2 | Classification & Group Summaries | 27 | 27% | | 3 | Resampling & Inference | 37 | 37% | | **Total** | | **100** | **100%** | : One point = one percent of the midterm grade. **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. --- # Exploring the Data *(36 points)* ## 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 {.unnumbered} **What to do** 1. Build a **grid of histograms**, one panel per numeric variable: `Age`, `Weight`, `Height`, and all ten measures in `wellness_vars` — thirteen panels. 2. Title each panel with its variable name. 3. 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. ```{r} #| fig-height: 7 # YOUR CODE HERE — Figure A ``` ![](F1.png){width=85%} ::: {.callout-note appearance="simple"} **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** a. Which distributions are roughly **symmetric**, and which have a long tail to the **right**? Name at least two of each. b. 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. c. 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 {.unnumbered} The ten measures are in different units, so you cannot put them on one axis as they stand. Standardise them first. **What to do** 1. 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. 2. Reshape to long format: one row per (measured variable name, z-score) pair. 3. **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. 4. **Order the measures by spread**, not alphabetically. 5. **Do not set axis limits.** If you do, some points will be dropped and R will warn you. Read the warning. ```{r} #| fig-height: 5 # YOUR CODE HERE — Figure B (boxplots) ``` ![](F2.png){width=85%} ::: {.callout-warning appearance="simple"} **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.]* --- ## A summary function, applied to every measure *(12 pts)* **What to do — part one: write the function** Write a function with exactly this signature: ```r wellness_summary <- function(df, var) { ... } ``` 1. Inside the function, extract the column named `var` from `df`. Refer to it programmatically — `df[[var]]` — not by typing the column name. 2. 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 | ```{r} wellness_summary <- function(df, var) { # YOUR CODE HERE } wellness_summary(df, "BMI") ``` **What to do — part two: apply it to all ten** 3. Call `wellness_summary()` on each of the ten measures in `wellness_vars`. 4. Assemble the results into a **numeric matrix with 7 rows and 10 columns** — one row per summary statistic, one column per measure. 5. 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. ```{r} # 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** a. 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. b. 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.]* --- ## The highest and lowest value of each measure *(8 pts)* **What to do** 1. For each of the ten measures in `wellness_vars`, find the respondent with the **highest** value and the respondent with the **lowest** value. 2. Collect the results in a data frame called `extremes` with exactly these five columns: `measure`, `max_id`, `max_value`, `min_id`, `min_value`. 3. 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`. ```{r} 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") ``` **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: a. Which measure, and what is the impossible value? b. How many respondents have it? (Count them — do not guess.) c. 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.]* --- # Classification & Group Summaries *(27 points)* ## A function that assigns a wellness profile *(14 pts)* **What to do** 1. Write a function with exactly this signature: ```r 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 | 2. **Apply your function to every respondent** in `df` to produce a character vector of length `nrow(df)`. A `for` loop over rows, or`sapply()`, will work. 3. 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. ```{r} 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" classify_wellness(145, 4) # expect "Elevated BP" (both rules -> first wins) classify_wellness(110, 4) # expect "Short Sleep" classify_wellness(110, 8) # expect "Neither" ``` ```{r} 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) ``` **What to report** Run these three checks and show the output: 1. every respondent received a label, and none is `NA`; 2. the **count of respondents satisfying both rules**; 3. How many respondents satisfy `sbp` ≥ 130 and `sleep` < 6? What are their wellness profile values? ```{r} # 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.]* --- ## Profile composition by age decade *(13 pts)* **What to do — part one: the overall counts** 1. Count the respondents in each profile, and compute each count as a **percentage** of the total. Display both in a small table. ```{r} 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** 2. For each age decade, find which profile is **most common**. 3. Build a data frame with columns `AgeDecade` and `dominant_profile`, and display it. 4. 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 | ```{r} 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") ``` **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:** 1. Build a table of counts: `table(df$AgeDecade, df$profile)`. 2. Convert to row-wise proportions with `prop.table(tab, margin = 1)`, so each decade sums to 1. 3. Convert to a data frame for plotting with `as.data.frame()`. You will get three columns; rename them `AgeDecade`, `Profile` and `Freq`. 4. Plot with `geom_col()`, multiplying `Freq` by 100 to get percentages. **Use this colour scheme** so everyone's figure is comparable: ```r scale_fill_manual(values = c("Elevated BP" = "#4C72B0", "Short Sleep" = "#DD8452", "Neither" = "#BFBFBF")) ``` A skeleton of the plotting call: ```r 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) ``` ```{r} # YOUR CODE HERE — build `comp`, then draw the stacked bar chart ``` ![](F3.png){width=85%} ::: {.callout-note appearance="simple"} **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** a. 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. b. Give the percentage in `"Elevated BP"` for the youngest decade and for the oldest. How large is the change? c. 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.]* --- # 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. ## 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** 1. Compute the observed statistic $T_{\text{obs}} = \bar{X}_{\text{inactive}} - \bar{X}_{\text{active}}$. 2. 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`. 3. Compute the two-sided p-value as $(r )/(B)$, where $r$ is the number of permuted differences with $|{\text{difference}}| \ge |T_{\text{obs}}|$. 4. Draw a **histogram** of `perm_diffs` with vertical lines at $T_{\text{obs}}$ and at $-T_{\text{obs}}$. 5. 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$. ```{r} 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]) ``` **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** a. State your conclusion. **Word it as a statement about the respondents in this sample**, not about Americans in general — see the dataset note above. b. 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.]* --- ## 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** 1. **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**. 2. **Is age related to the outcome?** Compute mean `BPSysAve` within each `AgeDecade`, and compute `cor(Age, BPSysAve)`. 3. **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. 4. **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. ```r w <- table(df$AgeDecade) / nrow(df) adjusted <- sum(w * within_decade_differences) ``` 5. Compare that adjusted figure with the unadjusted difference from 3.1, and compute what percentage of the original difference remains. ```{r} # Step 1: is age related to the grouping? # YOUR CODE HERE # Step 2: is age related to the outcome? # YOUR CODE HERE ``` ```{r} # Steps 3-4: within-decade differences, then the weighted average # YOUR CODE HERE ``` ```{r} # 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. A confounder must satisfy **two** conditions. State them, then show — using your own numbers from steps 1 and 2 — that age satisfies both here. b. How much of the unadjusted difference survives once age is held fixed? Give the two numbers and the percentage change. c. 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? d. 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.]* --- ## 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** 1. Compute the observed difference in means, $T_{\text{obs}} = \bar{X}_{\text{inactive}} - \bar{X}_{\text{active}}$. 2. 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`. 3. Report the **bootstrap standard error** — the standard deviation of `boot_diffs`. 4. Report a **95% percentile confidence interval** — the 2.5th and 97.5th percentiles of `boot_diffs`. 5. Draw a **histogram** of `boot_diffs` with vertical lines at $T_{\text{obs}}$ and at zero. 6. 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. ```{r} 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** a. Does the interval contain zero? State your conclusion, **worded as a statement about the respondents in this sample**. b. Compare your conclusion with the permutation p-value from 3.1. Do the two agree? c. 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? d. This interval is for the **unadjusted** difference. In light of 3.2, what should a reader be told alongside it? *[Write your answer here.]* --- # Declaration and References *(required)* ::: {.callout-important} ## 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 {.unnumbered} 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]* ::: {.callout-note appearance="simple"} 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 {.unnumbered} 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: 1. Stack Overflow, "How to order a factor by median before plotting." — accessed 2026-10-09. Used for the `reorder()` call in Figure B. 2. R documentation, `?quantile`. Consulted to check how ties are handled when computing the bootstrap percentile interval in 3.3. *[Your references here.]* ::: {.callout-tip appearance="simple"} **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 {.unnumbered} ```{r} sessionInfo() ```