diff --git a/NEWS.md b/NEWS.md index cedf22f..5f35620 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # primer.tutorials (development version) +* Fixed render failures in `dgm-logistic`'s Marginal Effects section under recent `marginaleffects`. `plot_predictions()` on the parsnip `logistic_reg()` fit either exhausts memory (`condition = `) or crashes (`by =`); and `predictions()` / `avg_predictions()` on the parsnip wrapper do not reliably echo the grid or `by` variable as a named column, so the `ggplot()` code fails with "object not found". The four plotting exercises (5, 10, 11, 12) now call `predictions()` / `avg_predictions()` on `extract_fit_engine(mod_logistic)` --- the underlying `glm`, whose `type = "response"` scale is the probability of leaving directly --- with an explicit `tibble()` / `expand_grid()` grid where a grid is needed. The printed exploratory chunks (1--4, 6--9) still use the tidymodels wrapper with `type = "prob"` and `group == "Yes"`. `plot_predictions()` is no longer used in the section. + +* **Split `05-recruits` into four per-virtue tutorials** --- `05-recruits-wisdom`, `05-recruits-justice`, `05-recruits-courage`, and `05-recruits-temperance` --- each with its own Introduction and Summary and its own YAML `id:`. Students carry one `analysis.qmd` in one `recruits` repo across the sequence. The `05-recruits-courage` tutorial absorbs the tidymodels content (Structure → Method → Fitting) from the retired `dgm-linear` tutorial, and `05-recruits-temperance` absorbs its marginaleffects content (Question → Units → Answer, the marginal/conditional distinction, `avg_predictions()`); `dgm-linear` is removed (its siblings `dgm-logistic`, `dgm-multinomial`, and `dgm-ordered` each re-teach the same framing and are unaffected). Knowledge drops in the Courage and Temperance tutorials are sourced from the DGM essays in `PPBDS/primer` (`guide/dgm.md`) or from the data. **Breaking change for student progress records:** completion data keyed on `05-recruits` or `dgm-linear` does not carry forward. + * Renamed `08-sps` to `08-seguro-popular` so the tutorial slug matches the program name students are studying (the dataset in `primer.data` is still named `sps`, but that internal tibble name was not a useful identifier for the tutorial). Tutorial directory, chapter file, YAML `id:`, and student repo name all change from `sps` to `seguro-popular`; the fitted-model variable changes from `fit_sps` to `fit_sp` (matching the program's "SP" abbreviation). **Breaking change for student progress records:** completion data keyed on `08-sps` does not carry forward. * **Renumbered all tutorials so that each tutorial's number matches its matching book chapter.** Specifically: discarded the `04-mechanics` tutorial (no longer in the curriculum); renamed `05-cardinal-virtues` → `04-cardinal-virtues`; and renumbered every example tutorial down by one (`06-recruits` → `05-recruits`, `07-trains` → `06-trains`, `08-colleges` → `07-colleges`, `09-sps` → `08-sps`, `10-smokes` → `09-smokes`, `11-shaming` → `10-shaming`). The matching book chapters in `book/` were renamed correspondingly. YAML `id:` fields and the `description:` "Tutorial #NN" lines were updated to match the new numbers. **Breaking change for student progress records:** completion data keyed on any of the old IDs (`04-mechanics`, `05-cardinal-virtues`, `06-recruits`, `07-trains`, `08-colleges`, `09-sps`, `10-smokes`, `11-shaming`) does not carry forward. The total curriculum is now 4 misc + 12 example = 16 tutorials, with the EMH split at `05–08` Easy, `09–12` Medium, `13–16` Hard. diff --git a/R/recruits.R b/R/recruits.R index 9295cf1..f012abd 100644 --- a/R/recruits.R +++ b/R/recruits.R @@ -1,11 +1,12 @@ #' Teaching cut of NHANES young-adult heights #' #' A 50-row sample drawn from the National Health and Nutrition Examination -#' Survey (NHANES, via [primer.data::nhanes]) for the `05-recruits` -#' tutorial. Restricted to ages 18--27 and to the columns relevant for the -#' tutorial's question. The split between sexes (40 male, 10 female) is +#' Survey (NHANES, via [primer.data::nhanes]) for the `05-recruits-*` +#' tutorials (`05-recruits-wisdom`, `-justice`, `-courage`, `-temperance`). +#' Restricted to ages 18--27 and to the columns relevant for the +#' tutorials' question. The split between sexes (40 male, 10 female) is #' deliberately uneven so the two group means have visibly different -#' standard errors --- a feature the tutorial's Temperance section asks +#' standard errors --- a feature the `05-recruits-temperance` tutorial asks #' students to notice and explain. #' #' @format A tibble with 50 rows and 3 variables: diff --git a/data-raw/recruits.R b/data-raw/recruits.R index dc46114..4c45bf2 100644 --- a/data-raw/recruits.R +++ b/data-raw/recruits.R @@ -1,9 +1,10 @@ # Build the `recruits` dataset shipped with the primer.tutorials package. # -# The 05-recruits tutorial uses this 50-row sample as a teaching cut of the -# NHANES height data. The deliberate mix --- 40 men and 10 women --- gives -# the two group means visibly different standard errors, which is one of the -# things the tutorial is meant to make students notice and explain. +# The 05-recruits-* tutorials (wisdom / justice / courage / temperance) use +# this 50-row sample as a teaching cut of the NHANES height data. The +# deliberate mix --- 40 men and 10 women --- gives the two group means +# visibly different standard errors, which is one of the things the +# 05-recruits-temperance tutorial is meant to make students notice and explain. # # Re-run this script (Rscript data-raw/recruits.R from the package root) # whenever you want to regenerate data/recruits.rda. usethis::use_data() diff --git a/inst/tutorials/03-rubin-causal-model/tutorial.Rmd b/inst/tutorials/03-rubin-causal-model/tutorial.Rmd index 7f934f9..a919e59 100644 --- a/inst/tutorials/03-rubin-causal-model/tutorial.Rmd +++ b/inst/tutorials/03-rubin-causal-model/tutorial.Rmd @@ -880,7 +880,7 @@ This tutorial walked through the Rubin Causal Model framework by applying the fo * **The fundamental problem of causal inference is that we only ever observe one potential outcome.** The other is the counterfactual, and we can never see it. * **The predictive/causal distinction is an analyst's commitment**, not a property of the data or the model. The same fit serves both questions; what changes is the language and the assumptions you are willing to defend. -The example tutorials that follow (`05-recruits`, `06-trains`, `07-colleges`, ...) walk through the same arc using real data and `tidymodels` code. The conceptual moves rehearsed here --- writing down the Preceptor Table, distinguishing potential outcomes from observable outcomes, defending unconfoundedness or acknowledging where it fails --- are exactly the moves you will use there. +The example tutorials that follow (`05-recruits-wisdom`, `06-trains`, `07-colleges`, ...) walk through the same arc using real data and `tidymodels` code. The conceptual moves rehearsed here --- writing down the Preceptor Table, distinguishing potential outcomes from observable outcomes, defending unconfoundedness or acknowledging where it fails --- are exactly the moves you will use there. The world is always more uncertain than our models would have us believe. diff --git a/inst/tutorials/04-cardinal-virtues/tutorial.Rmd b/inst/tutorials/04-cardinal-virtues/tutorial.Rmd index 3a2b16d..7453f8a 100644 --- a/inst/tutorials/04-cardinal-virtues/tutorial.Rmd +++ b/inst/tutorials/04-cardinal-virtues/tutorial.Rmd @@ -849,7 +849,7 @@ question_text(NULL, This tutorial walked through the four Cardinal Virtues --- *Wisdom*, *Justice*, *Courage*, and *Temperance* --- by applying them to two paired questions about high-school SAT scores: a predictive question (next year's average) and a causal one (the effect of a school-wide tutoring class). The exercises were conceptual, not coded. The takeaway is the framework itself: every data science problem proceeds through the same four-step arc, and the predictive/causal distinction is something the analyst commits to in Wisdom and defends in Justice --- not something the data dictates. -The example tutorials that follow (`05-recruits`, `06-trains`, `07-colleges`, ...) walk through the same arc using real data and `tidymodels` code. The hardest concepts in those tutorials --- writing down the Preceptor Table, defending the four assumptions, choosing the probability family --- are exactly the conceptual moves rehearsed here. +The example tutorials that follow (`05-recruits-wisdom`, `06-trains`, `07-colleges`, ...) walk through the same arc using real data and `tidymodels` code. The hardest concepts in those tutorials --- writing down the Preceptor Table, defending the four assumptions, choosing the probability family --- are exactly the conceptual moves rehearsed here. The world is always more uncertain than our models would have us believe. diff --git a/inst/tutorials/05-recruits-courage/tutorial.Rmd b/inst/tutorials/05-recruits-courage/tutorial.Rmd new file mode 100644 index 0000000..08330e4 --- /dev/null +++ b/inst/tutorials/05-recruits-courage/tutorial.Rmd @@ -0,0 +1,638 @@ +--- +title: "Recruits: Courage" +author: David Kane +tutorial: + id: 05-recruits-courage +output: + learnr::tutorial: + progressive: yes + allow_skip: yes +runtime: shiny_prerendered +description: "Tutorial #05 (Courage) for Preceptor's Primer" +--- + +```{r setup, include = FALSE} +library(learnr) +library(tutorial.helpers) +library(gt) + +library(tidyverse) +library(tidymodels) +library(primer.tutorials) +library(broom) + +knitr::opts_chunk$set(echo = FALSE) +knitr::opts_chunk$set(out.width = "90%") +options(tutorial.exercise.timelimit = 600, + tutorial.storage = "local") + +fit_recruits <- linear_reg() |> + set_engine("lm") |> + fit(height ~ sex, data = recruits) +``` + +```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} +``` + + +## Introduction +### + +This tutorial supports the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science: Using the Cardinal Virtues for Inference*](https://ppbds.github.io/primer/) by [David Kane](https://davidkane.info/). + +The world confronts us. Make decisions we must. + +Imagine that you are a logistics analyst at U.S. Marine Corps Supply, planning next year's bootcamp uniform order. The Quartermaster General has a big-picture goal --- every recruit kitted out at the right size on day one, within a fixed procurement budget --- and trusts you to figure out what to order and in what quantities. Many models would inform the order: distributions of height, weight, chest, foot length; the male-female mix in incoming recruits; the size-tail dynamics that drive how many extra-large uniforms to keep in inventory. This tutorial builds just one of them: the average height of male and female recruits. The number alone won't determine the order, but the height-by-sex split is one good input to the larger forecast. There are many decisions to make. + +### + +This is the third of four Recruits tutorials, one per Cardinal Virtue. **Wisdom** produced the question --- *What is the average height of male and female USMC recruits next year?* --- along with the Preceptor Table and a look at the data. **Justice** audited the assumptions behind the project and chose the functional form of the model: linear, with normally distributed errors. **Courage**, here, builds the data generating mechanism itself. + +Here is the Preceptor Table from Wisdom --- the smallest table such that, if every cell were filled in, the question would be easy to answer: + +```{r} +#| echo: false +#| results: asis +pre_tibble_recap <- tibble::tribble( + ~`Recruit`, ~`Height (cm)`, ~`Sex`, + "Maya Castillo", "150", "Female", + "Tyler Brooks", "170", "Male", + "...", "...", "...", + "Grace Sullivan", "160", "Female" +) + +pre_title_footnote_recap <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." +pre_units_footnote_recap <- "Each row is one USMC recruit expected to enlist during the coming year. Missing rows represent the other recruits not shown." +pre_outcome_footnote_recap <- "Height in centimeters, measured at enlistment." +pre_covariate_footnote_recap <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." + +pre_gt_html_recap <- gt::gt(pre_tibble_recap, id = "preceptor_tbl_courage") |> + gt::tab_header(title = "Preceptor Table") |> + gt::tab_spanner(label = "Unit", id = "unit_span", + columns = c(`Recruit`)) |> + gt::tab_spanner(label = "Outcome", id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pre_title_footnote_recap, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pre_units_footnote_recap, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pre_outcome_footnote_recap, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pre_covariate_footnote_recap, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #preceptor_tbl_courage .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pre_gt_html_recap, + "
\n", + "```\n", + sep = "" +) +``` + +The data you have is never the thing you care about. The data is a shadow cast by some underlying process --- the *data generating mechanism* --- and it is the process you actually want to understand. Our `recruits` sample is one such shadow. + +### Exercise 1 + +You should be connected to a repo named `recruits`, with an `analysis.qmd` you have been building through the Wisdom and Justice tutorials. If you are starting fresh here: create the repo and connect to it, create a Quarto document titled `"Recruits"` with yourself as author, save it as `analysis.qmd`, put `library(tidyverse)`, `library(tidymodels)`, `library(primer.tutorials)`, and `library(broom)` in a setup chunk with `#| message: false`, add + +``` +execute: + echo: false +``` + +to the YAML header, render, and open `analysis.html` with Live Server. + +In the R Terminal, run `show_file("analysis.qmd", chunk = "setup")`. CP/CR. + +```{r introduction-1} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +The `recruits` tibble is a 50-row teaching cut of NHANES: 40 male and 10 female young adults aged 18--27, with columns `height` (centimeters), `sex`, and `age`. The deliberately uneven 40/10 split will matter when we get to the precision of the two group means in Temperance. + +### Exercise 2 + +Add a code chunk to `analysis.qmd` that prints the `recruits` tibble (just `recruits` on its own line). Render. Copy the first few rows from the rendered HTML below. + +```{r introduction-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +Each row is one young adult, examined once. `height` is measured directly by an examiner; `sex` is recorded as Male or Female; `age` is in whole years. These 50 rows are the data --- the shadow. Courage's job is to write down a guess at the mechanism that cast it. + + +## Courage +### + +Justice gives us the Population Table and the abstract data generating mechanism. + +```{r} +#| echo: false +#| results: asis +pop_tibble_c <- tibble::tribble( + ~Source , ~`Young Adult` , ~`Year`, ~`Height (cm)`, ~`Sex` , + "..." , "..." , "..." , "..." , "..." , + "Data" , "Emily Tran" , "2002" , "162" , "Female" , + "Data" , "James Okafor" , "2005" , "178" , "Male" , + "Data" , "..." , "..." , "..." , "..." , + "Data" , "Tyler Brooks" , "2024" , "170" , "Male" , + "..." , "..." , "..." , "..." , "..." , + "Preceptor", "Maya Castillo" , "2026" , "150" , "Female" , + "Preceptor", "Tyler Brooks" , "2026" , "170" , "Male" , + "Preceptor", "..." , "..." , "..." , "..." , + "Preceptor", "Grace Sullivan" , "2026" , "160" , "Female" , + "..." , "..." , "..." , "..." , "..." +) + +pop_title_footnote_c <- "This table combines NHANES adult-height data from 1999 onward with the Preceptor Table's 2026 USMC recruits, drawn from the same broader population of US young adults aged 18--27." +pop_units_footnote_c <- "Each row is one young adult at one point in time. Data rows are NHANES survey participants between the ages of 18 and 27. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Tyler Brooks appears in both blocks." +pop_outcome_footnote_c <- "Height in centimeters. NHANES measures height directly in the Mobile Examination Center; USMC records height at enlistment." +pop_covariate_footnote_c <- "Sex, recorded as Male or Female in both data sources." + +pop_gt_html_c <- gt::gt(pop_tibble_c, id = "population_tbl_courage") |> + gt::tab_header(title = "Population Table") |> + gt::tab_spanner(label = "Unit/Time", id = "unit_span", + columns = c(`Young Adult`, `Year`)) |> + gt::tab_spanner(label = "Outcome" , id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Source`, `Young Adult`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Year`, `Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pop_title_footnote_c, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pop_units_footnote_c, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pop_outcome_footnote_c, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pop_covariate_footnote_c, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #population_tbl_courage .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pop_gt_html_c, + "
\n", + "```\n", + sep = "" +) +``` + +$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_n X_n + \epsilon$$ + +with $\epsilon \sim N(0, \sigma^2)$. + +### Exercise 1 + +In your own words, describe the components of the virtue of Courage for analyzing data. + +```{r courage-1} +question_text(NULL, + message = "Courage creates the data generating mechanism.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Building a model of the data generating mechanism takes three steps, and none of them has anything to do with software. + +First, **propose a structure**. What kind of outcome does the mechanism produce? The distribution of the outcome variable determines the structure of the DGM. A mechanism that produces numerical heights is a fundamentally different machine from one that produces pass/fail results. + +Second, **choose a method** for learning about the mechanism from data --- your philosophy of estimation. Will you seek the single version of the machine that best fits the data, or keep a whole collection of plausible versions? + +Third, **fit the DGM**. Your structure is really a family of machines --- regressions with every conceivable slope and intercept. Fitting means using the data to work out which settings are most plausible. + +**Structure → Method → Fitting.** That three-word chain is the discipline of the whole enterprise. + +### + +The abstract data generating mechanism shown in the preamble uses generic variables --- $Y$, $X_1$, $X_2$, and so on --- because its purpose is to describe the general mathematical structure. At the end of Justice we chose the functional family (linear with normal errors), but not the number or identity of the covariates we will ultimately include. Courage fills in the missing pieces. + +### Exercise 2 + +Because our outcome variable is continuous, add a code chunk to `analysis.qmd` that begins building the model with `linear_reg(engine = "lm")`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r courage-2-test} +#| echo: true +linear_reg(engine = "lm") +``` + +### + +`linear_reg()` is step one: it announces the **structure**, proposing that heights are centered on a value that shifts with the covariates, plus random noise. If our outcome recorded pass/fail instead, the same substantive question would demand a different machine. + +`engine = "lm"` is step two, the **method**: R's classical least-squares machinery, which returns the single best-fitting version of the mechanism. The structure and the method are kept separate on purpose, because the DGM (the truth out there) and our procedure for estimating it (the tool in here) are separate things. + +### Exercise 3 + +Extend the chunk in `analysis.qmd` to `linear_reg(engine = "lm") |> fit(height ~ sex, data = recruits)`. This is our first candidate model: outcome on the binary `sex`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-3} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r courage-3-test} +#| echo: true +linear_reg(engine = "lm") |> + fit(height ~ sex, data = recruits) +``` + +### + +`fit()` is step three, where the mechanism meets real data. You hand it a formula and the data, and the engine estimates the parameters --- the numerical settings of the proposed machine that best explain the shadow you have. + +A two-level character variable like `sex` is turned into a 0/1 "dummy" variable, renamed `sexMale` (the level for which the dummy is 1). Math doesn't speak text; dummies do the translation. + +### Exercise 4 + +Here is the parameter table for this candidate model: + +```{r} +#| echo: false +linear_reg() |> + set_engine("lm") |> + fit(height ~ sex, data = recruits) |> + tidy(conf.int = TRUE) |> + select(term, estimate, conf.low, conf.high) |> + mutate(across(where(is.numeric), \(x) round(x, 1))) +``` + +Write a sentence interpreting the `sexMale` estimate of 8.9. + +```{r courage-4} +question_text(NULL, + message = "When we compare two groups of recruits differing only in sex, the male group has an expected height about 8.9 centimeters taller than the female group. The intercept of 165.9 is the expected height of a female recruit; adding `sexMale` (= 1 for male) gives the expected male height of about 174.8 cm.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +Notice the language: *"when we compare two groups,"* not *"raising sex by one unit increases height."* For a binary covariate the comparison framing is the only honest reading. The 95% confidence interval for `sexMale` is [4.3, 13.4] centimeters; zero is well outside that interval, so the model says the *expected* height of male recruits is taller than the expected height of female recruits. That does not mean every male is taller than every female --- the two distributions overlap. + +### Exercise 5 + +Change the formula to `fit(height ~ age, data = recruits)`. This refits with `age`, a continuous covariate, instead of `sex`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-5} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r courage-5-test} +#| echo: true +linear_reg(engine = "lm") |> + fit(height ~ age, data = recruits) +``` + +### + +Continuous covariates produce a single slope coefficient. The reading is a per-unit comparison: how much do two groups, differing by one unit of the covariate, differ in expected outcome? + +### Exercise 6 + +Here is the parameter table for the second candidate model: + +```{r} +#| echo: false +linear_reg() |> + set_engine("lm") |> + fit(height ~ age, data = recruits) |> + tidy(conf.int = TRUE) |> + select(term, estimate, conf.low, conf.high) |> + mutate(across(where(is.numeric), \(x) round(x, 1))) +``` + +Write a sentence interpreting the slope on `age`. + +```{r courage-6} +question_text(NULL, + message = "The slope on `age` is essentially zero (0.0 cm per year, with a 95% confidence interval that crosses zero in both directions). When we compare two groups of recruits differing by one year of age (within our 18--27 window), their expected heights are indistinguishable. Adult height has stopped changing by the late teens; this fit confirms that.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +There is no single "correct" DGM for a given question. Which right-hand-side variables to include is a judgment call --- one we make, and could remake. We tried `age`, and its wide confidence interval, straddling zero in both directions, is exactly what *"this covariate doesn't help"* looks like in a parameter table. Our final model is therefore `height ~ sex` --- the candidate from Exercise 4 --- which is the model the rest of the tutorial uses. + +### Exercise 7 + +We've settled on the final fit. In `analysis.qmd`, add a new code chunk that assigns the result of `linear_reg() |> set_engine("lm") |> fit(height ~ sex, data = recruits)` to `fit_recruits`, then on the next line prints `fit_recruits`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-7} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r courage-7-test} +#| echo: true +fit_recruits <- linear_reg() |> + set_engine("lm") |> + fit(height ~ sex, data = recruits) +fit_recruits +``` + +### + +Only once the parameters are estimated do you have a usable DGM --- one you can use to predict new observations, or to ask what would happen if you changed one part of the mechanism. `fit_recruits` is that object: the fitted machine every later question will be addressed to. + +### Exercise 8 + +Re-fitting the model on every render is wasteful --- the fit's the same each time. `fit_recruits` earns its own cached create chunk. The next three exercises set that up: clean the chunk, turn caching on, add the cache directory to `.gitignore`. + +First, clean up. In your chunk from Exercise 7, delete the `fit_recruits` line at the bottom, leaving just the `fit_recruits <- linear_reg() ...` assignment. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-8} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 8) +``` + +### + +
fit_recruits <- linear_reg() |>
+  set_engine("lm") |>
+  fit(height ~ sex, data = recruits)
+
+ +### + +A chunk ready to cache should contain only the artifact it builds --- nothing else. + +### Exercise 9 + +Add `#| cache: true` as the first line inside the chunk. Render the QMD. Rendering creates an `analysis_cache/` directory next to your `analysis.qmd`, which holds the saved fit so Quarto can reload it on subsequent renders rather than re-fitting every time. To confirm the directory is there, run `ls` in the **bash Terminal**. CP/CR. + +```{r courage-9} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +
$ ls
+analysis.qmd  analysis_cache  analysis_files  analysis.html
+
+ +### + +For a 50-row linear regression the savings are small; for the larger model fits in later tutorials, caching can be the difference between a one-second and a one-minute render. + +### Exercise 10 + +The Source Control change count just jumped, because `analysis_cache/` and everything inside it now look like untracked work to git. Cached files should never go to GitHub --- they're machine-specific and regenerated on every render. + +Add `analysis_cache` to your `.gitignore` on its own line. In the R Terminal, run `show_file(".gitignore")`. CP/CR. + +```{r courage-10} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +
analysis_files
+analysis_cache
+
+ +### + +The Source Control change count should drop back to whatever it was before the render. Cached objects live next to the QMD, not in version control. + +With the parameters estimated, we can write the fitted model as a concrete formula: + +$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ + +with residuals drawn from $N(0, 41.0)$ --- a residual standard deviation of about 6.4 cm. + +Three differences from the abstract form in the preamble. First, the parameters are replaced by their estimates. Second, the error term is gone --- this version of the DGM answers an expected-value question, so the randomness sits to the side in the residuals line until a question needs it. Third, the left-hand side has a hat, because hats mark estimated values. + +**This is our data generating mechanism.** A data generating mechanism is a formal, mathematical guess about the process that cast the shadow --- something we can write down and implement with computer code. The randomness built into it lives in the residuals line. + +### Exercise 11 + +Start a new working chunk below the cached `fit_recruits` chunk that runs `tidy(fit_recruits, conf.int = TRUE)`. This returns 95% intervals for all the parameters in our model. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r courage-11} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r courage-11-test} +#| echo: true +tidy(fit_recruits, conf.int = TRUE) +``` + +### + +Here is a nicer-looking version of that parameter table: + +```{r} +#| echo: false +tidy(fit_recruits, conf.int = TRUE) |> + select(term, estimate, conf.low, conf.high) |> + mutate(across(where(is.numeric), \(x) round(x, 1))) |> + knitr::kable( + caption = "Estimates and 95% Confidence Intervals for Parameters. Source: 50-row sample from NHANES via primer.tutorials::recruits" + ) +``` + +Rounding to one decimal place keeps the displayed values consistent with the fitted formula above and is fine enough that the estimate and its confidence-interval bounds stay visibly distinct in every row --- you need rounding finer than the standard error to keep adjacent values apart. At the very least, a table like this should carry a title and a caption naming the data source. + +### + +We don't really care about coefficients. Nobody ever asked, *"What are the parameters of a linear regression relating height to sex?"* The fitted DGM is not the answer to anything --- it is the machine we built so that we can answer questions like *"What height should we expect for a male recruit?"* The Temperance tutorial runs that machine. + +### Exercise 12 + +Add a sentence to your project summary explaining the structure of the model. Something like: "We model XX [concept of outcome, not variable name], [values of XX], as a [linear/logistic/multinomial/ordinal] function of XX [and maybe other covariates]." + +Recall the beginning of our version of the summary: + +> People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. + +```{r courage-12} +question_text(NULL, + message = "We model height as a normally distributed variable which is a linear function of sex.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Read our answer. Do not copy/paste exactly. Add your two sentences to the summary paragraph. Render the QMD. + +### Exercise 13 + +Commit and push, with a message like "Fit final linear regression DGM for recruits." In the bash Terminal, run: + +``` +git log -1 +``` + +CP/CR. + +```{r courage-13} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +`fit_recruits` is saved, cached, and committed. It is not the answer to anything by itself --- it is the machine we built so we can answer questions, which is exactly what the Temperance tutorial does with it. + + +## Summary +### + +This tutorial covered the **Courage** stage of the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science*](https://ppbds.github.io/primer/). + +Building a model of the data generating mechanism takes three steps: **propose a structure** (the distribution of the outcome variable determines it --- a continuous height calls for `linear_reg()`), **choose a method** (`set_engine("lm")` for a classical least-squares fit), and **fit the DGM** (`fit(height ~ sex, data = recruits)`, which estimates the parameters). Structure → Method → Fitting. + +We tried `sex` and `age` as candidate covariates, kept `sex` because its confidence interval cleared zero and dropped `age` because its did not, and saved the final fit as `fit_recruits`: + +$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ + +with residuals drawn from $N(0, 41.0)$. In words: we model height as a normally distributed variable which is a linear function of sex. + +The fitted DGM is not the answer to anything by itself. It is the machine we built so that we can answer questions --- expected heights, and the uncertainty around them --- which the **Temperance** tutorial takes up next. + +The world is always more uncertain than our models would have us believe. + +```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} +``` diff --git a/inst/tutorials/05-recruits-justice/tutorial.Rmd b/inst/tutorials/05-recruits-justice/tutorial.Rmd new file mode 100644 index 0000000..54f92fe --- /dev/null +++ b/inst/tutorials/05-recruits-justice/tutorial.Rmd @@ -0,0 +1,628 @@ +--- +title: "Recruits: Justice" +author: David Kane +tutorial: + id: 05-recruits-justice +output: + learnr::tutorial: + progressive: yes + allow_skip: yes +runtime: shiny_prerendered +description: "Tutorial #05 (Justice) for Preceptor's Primer" +--- + +```{r setup, include = FALSE} +library(learnr) +library(tutorial.helpers) +library(gt) + +library(tidyverse) +library(primer.tutorials) + +knitr::opts_chunk$set(echo = FALSE) +knitr::opts_chunk$set(out.width = "90%") +options(tutorial.exercise.timelimit = 600, + tutorial.storage = "local") +``` + +```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} +``` + + +## Introduction +### + +This tutorial supports the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science: Using the Cardinal Virtues for Inference*](https://ppbds.github.io/primer/) by [David Kane](https://davidkane.info/). + +The world confronts us. Make decisions we must. + +Imagine that you are a logistics analyst at U.S. Marine Corps Supply, planning next year's bootcamp uniform order. The Quartermaster General has a big-picture goal --- every recruit kitted out at the right size on day one, within a fixed procurement budget --- and trusts you to figure out what to order and in what quantities. Many models would inform the order: distributions of height, weight, chest, foot length; the male-female mix in incoming recruits; the size-tail dynamics that drive how many extra-large uniforms to keep in inventory. This tutorial builds just one of them: the average height of male and female recruits. The number alone won't determine the order, but the height-by-sex split is one good input to the larger forecast. There are many decisions to make. + +### + +This is the second of four Recruits tutorials, one per Cardinal Virtue. **Wisdom** produced the question --- *What is the average height of male and female USMC recruits next year?* --- along with the Preceptor Table and a first look at the `recruits` data. **Justice**, here, audits the assumptions that let us treat the data and the Preceptor Table as two views of one population, and then chooses the functional form of the model. + +Here is the Preceptor Table from Wisdom --- the smallest table such that, if every cell were filled in, the question would be easy to answer: + +```{r} +#| echo: false +#| results: asis +pre_tibble_recap <- tibble::tribble( + ~`Recruit`, ~`Height (cm)`, ~`Sex`, + "Maya Castillo", "150", "Female", + "Tyler Brooks", "170", "Male", + "...", "...", "...", + "Grace Sullivan", "160", "Female" +) + +pre_title_footnote_recap <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." +pre_units_footnote_recap <- "Each row is one USMC recruit expected to enlist during the coming year. Missing rows represent the other recruits not shown." +pre_outcome_footnote_recap <- "Height in centimeters, measured at enlistment." +pre_covariate_footnote_recap <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." + +pre_gt_html_recap <- gt::gt(pre_tibble_recap, id = "preceptor_tbl_recap") |> + gt::tab_header(title = "Preceptor Table") |> + gt::tab_spanner(label = "Unit", id = "unit_span", + columns = c(`Recruit`)) |> + gt::tab_spanner(label = "Outcome", id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pre_title_footnote_recap, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pre_units_footnote_recap, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pre_outcome_footnote_recap, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pre_covariate_footnote_recap, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #preceptor_tbl_recap .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pre_gt_html_recap, + "
\n", + "```\n", + sep = "" +) +``` + +### Exercise 1 + +You should be connected to a repo named `recruits`, with an `analysis.qmd` you started in the Wisdom tutorial. If you are starting fresh here: create the repo and connect to it, create a Quarto document titled `"Recruits"` with yourself as author, save it as `analysis.qmd`, put `library(tidyverse)` and `library(primer.tutorials)` in a setup chunk with `#| message: false`, add + +``` +execute: + echo: false +``` + +to the YAML header, render, and open `analysis.html` with Live Server. + +In the R Terminal, run `show_file("analysis.qmd", chunk = "setup")`. CP/CR. + +```{r introduction-1} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +The `recruits` tibble is a 50-row teaching cut of NHANES: 40 male and 10 female young adults aged 18--27, with columns `height` (centimeters), `sex`, and `age`. + +### Exercise 2 + +Add a code chunk to `analysis.qmd` that runs `recruits |> count(sex)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r introduction-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r introduction-2-test} +#| echo: true +recruits |> count(sex) +``` + +### + +40 male rows and 10 female. NHANES intentionally oversamples some demographic groups so each subgroup can be studied on its own, and our `recruits` sample amplifies that: we drew 40 men and 10 women rather than matching their representation in NHANES. Both are representativeness concerns, which the exercises below name. + + +## Justice +### + +Wisdom gives us the Preceptor Table and the data. + +```{r} +#| echo: false +#| results: asis +pre_tibble_j <- tibble::tribble( + ~`Recruit`, ~`Height (cm)`, ~`Sex`, + "Maya Castillo", "150", "Female", + "Tyler Brooks", "170", "Male", + "...", "...", "...", + "Grace Sullivan", "160", "Female" +) + +pre_title_footnote_j <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." +pre_units_footnote_j <- "Each row is one USMC recruit expected to enlist during the coming year. Missing rows represent the other recruits not shown." +pre_outcome_footnote_j <- "Height in centimeters, measured at enlistment." +pre_covariate_footnote_j <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." + +pre_gt_html_j <- gt::gt(pre_tibble_j, id = "preceptor_tbl_justice") |> + gt::tab_header(title = "Preceptor Table") |> + gt::tab_spanner(label = "Unit", id = "unit_span", + columns = c(`Recruit`)) |> + gt::tab_spanner(label = "Outcome", id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pre_title_footnote_j, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pre_units_footnote_j, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pre_outcome_footnote_j, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pre_covariate_footnote_j, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #preceptor_tbl_justice .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pre_gt_html_j, + "
\n", + "```\n", + sep = "" +) +``` + + +```{r} +#| echo: false +#| results: asis +data_tibble <- tibble::tribble( + ~`Respondent`, ~`Age`, ~`Height (cm)`, ~`Sex`, + "Emily Tran", "21", "162", "Female", + "James Okafor", "25", "178", "Male", + "Rachel Goldberg", "19", "159", "Female", + "...", "...", "...", "...", + "Tyler Brooks", "23", "170", "Male" +) + +data_title_footnote <- "National Health and Nutrition Examination Survey (NHANES), Centers for Disease Control and Prevention. Continuous survey data, 1999--present. We restrict to respondents aged 18 through 27 with non-missing height and sex --- about 1,300 rows. The original survey used a complex, stratified, multistage sampling design; we ignore the survey weights." +data_units_footnote <- "Each row is one NHANES participant examined between 1999 and the most recent release cycle. NHANES public files are anonymized, so the example rows use invented but plausible names." +data_outcome_footnote <- "Height in centimeters, measured by a trained examiner in NHANES's Mobile Examination Center using a fixed stadiometer. Each participant is measured once per survey visit with shoes removed." +data_covariate_footnote <- "Sex, recorded by the NHANES interviewer as 'Male' or 'Female'. NHANES does not separately record gender identity; the two-category coding dates from survey waves going back decades." + +data_gt_html <- gt::gt(data_tibble, id = "data_tbl") |> + gt::tab_header(title = "Data: NHANES Young Adults, Ages 18–27") |> + gt::tab_spanner(label = "Unit", id = "unit_span", + columns = c(`Respondent`, `Age`)) |> + gt::tab_spanner(label = "Outcome", id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Respondent`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Age`, `Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = data_title_footnote, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = data_units_footnote, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = data_outcome_footnote, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = data_covariate_footnote, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #data_tbl .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + data_gt_html, + "
\n", + "```\n", + sep = "" +) +``` + +### Exercise 1 + +In your own words, describe the components of Justice when working on a data science problem. + +```{r justice-1} +question_text(NULL, + message = "Justice reviews the Population Table and selects the formula for the data generating mechanism.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Justice is about concerns that you (or your critics) might have, reasons why the model you create might not work as well as you hope. The exercises that follow walk through several such concerns by name. Justice also picks the formula for the data generating mechanism --- different outcome variables call for different formulas, a piece that gets more technical names in later tutorials. + +### Exercise 2 + +In your own words, define "validity" as we use the term. + +```{r justice-2} +question_text(NULL, + message = "Validity is the consistency, or lack thereof, in the columns of the data set and the corresponding columns in the Preceptor Table.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Validity is always about the columns in the Preceptor Table and the data. Just because columns from these two different tables have the same name does not mean that they are the same thing. + +### Exercise 3 + +Provide one reason why the assumption of validity might not hold for the outcome variable `height` or for one of the covariates. Use the words "column" or "columns" in your answer. + +```{r justice-3} +question_text(NULL, + message = "The `height` column in the data is recorded by NHANES's trained examiners using a fixed stadiometer, with shoes removed. The `height` column in the Preceptor Table is recorded at USMC enlistment with whatever instrument the recruiting office has on hand --- possibly with shoes on, possibly under different instructions. Even though both columns are labeled 'Height (cm),' the two measurement procedures are not the same, so the columns may not represent the same underlying quantity.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +In order to consider the Preceptor Table and the data to be drawn from the same population, the columns from one must have a *valid correspondence* with the columns in the other. Validity, if true (or at least reasonable), allows us to construct the Population Table, which is the first step in Justice. + +### Exercise 4 + +In your own words, define a Population Table. + +```{r justice-4} +question_text(NULL, + message = "The Population Table includes a row for each unit/time combination in the underlying population from which both the Preceptor Table and the data are drawn.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +The Population Table is almost always much bigger than the combination of the Preceptor Table and the data, because if we can really assume that both are part of the same population, then that population must cover a broad universe of time and units. + +### Exercise 5 + +Specify the unit/time combinations which define each row in this Population Table. + +```{r justice-5} +question_text(NULL, + message = "Each row is one young adult (aged 18--27) in a given year. Data rows are NHANES participants examined between 1999 and the most recent release. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Other rows cover all the young adults alive during the overlapping span who appear in neither table.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +```{r} +#| echo: false +#| results: asis +pop_tibble <- tibble::tribble( + ~Source , ~`Young Adult` , ~`Year`, ~`Height (cm)`, ~`Sex` , + "..." , "..." , "..." , "..." , "..." , + "Data" , "Emily Tran" , "2002" , "162" , "Female" , + "Data" , "James Okafor" , "2005" , "178" , "Male" , + "Data" , "..." , "..." , "..." , "..." , + "Data" , "Tyler Brooks" , "2024" , "170" , "Male" , + "..." , "..." , "..." , "..." , "..." , + "Preceptor", "Maya Castillo" , "2026" , "150" , "Female" , + "Preceptor", "Tyler Brooks" , "2026" , "170" , "Male" , + "Preceptor", "..." , "..." , "..." , "..." , + "Preceptor", "Grace Sullivan" , "2026" , "160" , "Female" , + "..." , "..." , "..." , "..." , "..." +) + +pop_title_footnote <- "This table combines NHANES adult-height data from 1999 onward with the Preceptor Table's 2026 USMC recruits, drawn from the same broader population of US young adults aged 18--27." +pop_units_footnote <- "Each row is one young adult at one point in time. Data rows are NHANES survey participants between the ages of 18 and 27. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Tyler Brooks appears in both blocks --- surveyed by NHANES in 2024 at age 23, enlisting in 2026 at 25. Same person, different times, different rows." +pop_outcome_footnote <- "Height in centimeters. NHANES measures height directly in the Mobile Examination Center; USMC records height at enlistment. The two measurement procedures are close enough that we treat them as the same column, subject to the validity concern raised in Exercise 3." +pop_covariate_footnote <- "Sex, recorded as Male or Female in both data sources." + +pop_gt_html <- gt::gt(pop_tibble, id = "population_tbl") |> + gt::tab_header(title = "Population Table") |> + gt::tab_spanner(label = "Unit/Time", id = "unit_span", + columns = c(`Young Adult`, `Year`)) |> + gt::tab_spanner(label = "Outcome" , id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Source`, `Young Adult`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Year`, `Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pop_title_footnote, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pop_units_footnote, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pop_outcome_footnote, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pop_covariate_footnote, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #population_tbl .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pop_gt_html, + "
\n", + "```\n", + sep = "" +) +``` + +### Exercise 6 + +In your own words, define the assumption of "stability" when employed in the context of data science. + +```{r justice-6} +question_text(NULL, + message = "Stability means that the relationship between the columns in the Population Table is the same for three categories of rows: the data, the Preceptor Table, and the larger population from which both are drawn.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +The most natural way for stability to fail is *time*. Our data was gathered in earlier NHANES cycles, but our recruits will enlist next year. Stability holds only if the world that produced the data behaves the same as the world the recruits enter. That's not guaranteed here: NHANES spans 25 years, and compared to next year's recruit cohort, that's a long enough window for the relationship between sex and adult height to drift. +We will sharpen this idea --- the *parameters* of the model, not the levels of any single variable --- in the next exercise. + +### Exercise 7 + +Provide one reason why the assumption of stability might not be true in this case. + +```{r justice-7} +question_text(NULL, + message = "The *intercept* of the height-on-sex regression --- the expected height of a female young adult --- has likely drifted upward over the 25-year span of NHANES, while the *slope* (the male-female height gap, the `sexMale` coefficient) has stayed roughly constant. By 2026 the underlying $\\beta_0$ for our recruits is probably a centimeter or two higher than the average across NHANES's reporting years. That is a *parameter* shift; if it is large enough, stability is violated. The fact that average heights in the population have changed is not by itself the issue --- a distribution change in `height` is something a stable model can absorb. What stability requires is that $\\beta_0$ and $\\beta_1$ stay put.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +A common confusion is to point at any change between two eras and call it a stability violation. It isn't. The distribution of `height` in the population may have shifted upward; the spread of heights may have narrowed; the mix of ethnicities in the recruiting pool may differ from NHANES participants. None of those, on their own, is a stability violation. What stability requires is that the *parameters* of our model --- $\beta_0$ and $\beta_1$ --- be the same in the data, in the Preceptor Table, and in the broader population. Distribution shifts are everywhere; parameter shifts are what hurt us. + +### Exercise 8 + +We use our data to make inferences about the overall population. We use information about the population to make inferences about the Preceptor Table: Data → Population → Preceptor Table. In your own words, define the assumption of "representativeness." + +```{r justice-8} +question_text(NULL, + message = "Representativeness, or the lack thereof, concerns two relationships among the rows in the Population Table. The first is between the data and the other rows. The second is between the other rows and the Preceptor Table.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Ideally, we would like both the Preceptor Table and our data to be random samples from the population. Sadly, this is almost never the case. + +### Exercise 9 + +We do not use the data directly to estimate missing values in the Preceptor Table. Instead, we use the data to learn about the overall population. Provide one reason, involving the relationship between the data and the population, why the assumption of representativeness might not be true in this case. + +```{r justice-9} +question_text(NULL, + message = "NHANES is a voluntary survey. People who agree to participate (and who show up for the Mobile Examination Center visit) are not a random sample of US young adults --- they skew toward people with time, transportation, and some trust in government research. On top of that, the `recruits` tibble we are fitting on is itself a non-random subsample of NHANES: 40 male rows and 10 female rows, deliberately unbalanced. Both layers --- the survey participation and our 40/10 choice --- mean the data we model is not a random sample of US young adults.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +When representativeness is violated, our parameter estimates *might* be biased --- not "will be," because by pure luck a non-representative sample can still produce the right number, but we have no principled reason to expect that to happen. There is a *second* representativeness relationship --- between the population and the Preceptor Table --- that later tutorials will exercise; for now, the data ↔ population link is what we focus on. + +### Exercise 10 + +A statistical model consists of two parts: the probability family and the link function. The probability family is the probability distribution that generates the randomness in our data. Our outcome variable is `height` --- continuous. What probability family should we use? + +```{r justice-10} +question_text(NULL, + message = "Normal: $Y \\sim N(\\mu, \\sigma^2)$. The Normal distribution is the standard choice for continuous outcomes.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 4) +``` + +### + +While you're here, also add `library(tidymodels)` to the setup chunk in `analysis.qmd` so the rest of Courage can call `linear_reg()` and `fit()`. + +### Exercise 11 + +The link function is the mathematical formula that links our data to the unknown parameters. What is the link function for a continuous outcome? + +```{r justice-11} +question_text(NULL, + message = "Linear: $\\mu = \\beta_0 + \\beta_1 X_1 + \\ldots + \\beta_n X_n$. Our model will be a linear regression.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 4) +``` + +### + +While you're here, also add `library(broom)` to the setup chunk in `analysis.qmd` so we can call `tidy()` on the fitted model in Courage. + +### Exercise 12 + +Write one sentence highlighting a potential weakness in your model. Derive it from possible problems with the assumptions above. We will add this to our summary paragraph. So far our version of the summary paragraph looks like this: + +> People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year. + +Your version will be somewhat different. + +```{r justice-12} +question_text(NULL, + message = "One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Add a weakness sentence to the summary paragraph in your QMD. You can modify your paragraph, but don't copy/paste our answer exactly. Render the QMD, then commit/push. + + +## Summary +### + +This tutorial covered the **Justice** stage of the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science*](https://ppbds.github.io/primer/). + +Justice audits the assumptions that let us treat the data and the Preceptor Table as two views of one population: + +- **Validity** --- do the `height` and `sex` columns mean the same thing in the data (NHANES examiners, fixed stadiometer) as in the Preceptor Table (USMC enlistment)? +- **Stability** --- are the *parameters* of the height-on-sex regression the same across the 25-year span of NHANES and next year's recruit cohort? +- **Representativeness** --- is the data a fair draw from the population, given that NHANES is voluntary and our 40/10 male-female split is deliberately unbalanced? + +Justice also chose the model's functional form: because `height` is continuous, a Normal probability family with a linear link --- an ordinary linear regression, written abstractly as + +$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_n X_n + \epsilon$$ + +with $\epsilon \sim N(0, \sigma^2)$. The **Courage** tutorial fills in the covariates and estimates the parameters. + +The world is always more uncertain than our models would have us believe. + +```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} +``` diff --git a/inst/tutorials/05-recruits-temperance/tutorial.Rmd b/inst/tutorials/05-recruits-temperance/tutorial.Rmd new file mode 100644 index 0000000..7eca66e --- /dev/null +++ b/inst/tutorials/05-recruits-temperance/tutorial.Rmd @@ -0,0 +1,536 @@ +--- +title: "Recruits: Temperance" +author: David Kane +tutorial: + id: 05-recruits-temperance +output: + learnr::tutorial: + progressive: yes + allow_skip: yes +runtime: shiny_prerendered +description: "Tutorial #05 (Temperance) for Preceptor's Primer" +--- + +```{r setup, include = FALSE} +library(learnr) +library(tutorial.helpers) + +library(tidyverse) +library(tidymodels) +library(primer.tutorials) +library(broom) +library(marginaleffects) + +knitr::opts_chunk$set(echo = FALSE) +knitr::opts_chunk$set(out.width = "90%") +options(tutorial.exercise.timelimit = 600, + tutorial.storage = "local") + +fit_recruits <- linear_reg() |> + set_engine("lm") |> + fit(height ~ sex, data = recruits) +``` + +```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} +``` + + +## Introduction +### + +This tutorial supports the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science: Using the Cardinal Virtues for Inference*](https://ppbds.github.io/primer/) by [David Kane](https://davidkane.info/). + +The world confronts us. Make decisions we must. + +Imagine that you are a logistics analyst at U.S. Marine Corps Supply, planning next year's bootcamp uniform order. The Quartermaster General has a big-picture goal --- every recruit kitted out at the right size on day one, within a fixed procurement budget --- and trusts you to figure out what to order and in what quantities. Many models would inform the order: distributions of height, weight, chest, foot length; the male-female mix in incoming recruits; the size-tail dynamics that drive how many extra-large uniforms to keep in inventory. This tutorial builds just one of them: the average height of male and female recruits. The number alone won't determine the order, but the height-by-sex split is one good input to the larger forecast. There are many decisions to make. + +### + +This is the last of four Recruits tutorials, one per Cardinal Virtue. **Wisdom** produced the question and the Preceptor Table; **Justice** audited the assumptions and chose a linear model with normal errors; **Courage** fit the data generating mechanism and saved it as `fit_recruits`: + +$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ + +with residuals drawn from $N(0, 41.0)$ --- a residual standard deviation of about 6.4 cm. In words: we model height as a normally distributed variable which is a linear function of sex. + +Here is the summary paragraph so far, built up across the first three tutorials: + +> People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. We model height as a normally distributed variable which is a linear function of sex. + +**Temperance** interprets that mechanism and uses it --- with humility --- to answer the question. + +### Exercise 1 + +You should be connected to a repo named `recruits`, with an `analysis.qmd` you have been building through the Wisdom, Justice, and Courage tutorials. Its setup chunk should have `library(tidyverse)`, `library(tidymodels)`, `library(primer.tutorials)`, and `library(broom)`, and it should contain a cached chunk that assigns `fit_recruits`. + +In the R Terminal, run `show_file("analysis.qmd", chunk = "fit_recruits")` (or, if your chunk is unlabeled, `show_file("analysis.qmd")` and find the `fit_recruits` chunk). CP/CR. + +```{r introduction-1} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +
#| cache: true
+fit_recruits <- linear_reg() |>
+  set_engine("lm") |>
+  fit(height ~ sex, data = recruits)
+
+ +### Exercise 2 + +Confirm the fitted mechanism is what Courage produced. Add a chunk that runs `tidy(fit_recruits, conf.int = TRUE)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r introduction-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r introduction-2-test} +#| echo: true +tidy(fit_recruits, conf.int = TRUE) +``` + +### + +The intercept is about 165.9 cm (the expected height of a female recruit) and the `sexMale` coefficient is about 8.9 cm, with a 95% interval of roughly [4.3, 13.4]. These are the parameter values Temperance now puts to work. + + +## Temperance +### + +Courage provides the data generating mechanism. We can express the DGM in four ways: + +**In words:** We model height as a normally distributed variable which is a linear function of sex. + +**In R code:** + +```r +fit_recruits <- linear_reg() |> + set_engine("lm") |> + fit(height ~ sex, data = recruits) +``` + +**In a parameter table:** + +```{r} +#| echo: false +fit_recruits |> + tidy(conf.int = TRUE) |> + select(term, estimate, conf.low, conf.high) |> + mutate(across(where(is.numeric), \(x) round(x, 1))) +``` + +**As a mathematical formula:** + +$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ + +with residuals drawn from $N(0, 41.0)$ --- a residual standard deviation of about 6.4 cm. The hatted equation is the expected-value mode of the DGM --- plug in a sex, get an expected height --- while the residuals line keeps the draw mode's randomness at hand for questions about individual recruits. + +### Exercise 1 + +In your own words, describe the use of Temperance in data science. + +```{r temperance-1} +question_text(NULL, + message = "Temperance interprets the data generating mechanism and then uses it to answer, with the help of graphics, the question(s) with which we began. Humility reminds us that this answer is always false.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +The fitted DGM is not the answer to anything. Nobody ever asked, *"What are the parameters of a linear regression relating height to sex?"* The questions people actually ask sound like: *What height should we expect for a male recruit? How much taller, on average, are male recruits than female recruits?* + +Answering a question with a DGM takes three conceptual moves. First, **pose the question as a quantity** the machine can produce --- a *prediction* (the expected outcome for a specified unit) or a *comparison* (run the machine twice and report the difference). Second, **specify the units** --- one imaginary recruit, every recruit in the sample, a grid of hypothetical recruits. Third, **run the machine and summarize**, uncertainty included. **Question → Units → Answer.** The tool for these moves is the **[marginaleffects](https://marginaleffects.com/)** package. + +### Exercise 2 + +In the end, we don't really care about parameters, much less how to interpret them. Parameters are *imaginary*, like unicorns. We care about answers to our questions. *In the modern world, all parameters are nuisance parameters.* + +If you do not already have **marginaleffects** installed, run `install.packages("marginaleffects")` in the R Terminal. Then add `library(marginaleffects)` to the setup chunk in `analysis.qmd` and render. In the R Terminal, run `show_file("analysis.qmd", chunk = "setup")`. CP/CR. + +```{r temperance-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 8) +``` + +### + +

+#| message: false
+library(tidyverse)
+library(tidymodels)
+library(primer.tutorials)
+library(broom)
+library(marginaleffects)
+
+ +### + +Be modest in the claims you make. When you state a conclusion, ask whether a different, equally plausible mechanism could have produced the same data. If it could, your conclusion is weaker than it looks. Yet decisions made with a flawed model are almost always better than decisions made without one. + +### Exercise 3 + +What is the specific question we are trying to answer? + +```{r temperance-3} +question_text(NULL, + message = "What is the expected height of male and female USMC recruits next year?", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +Vague curiosity --- *"does sex matter for height?"* --- has to become a number the machine can produce. Our question asks for two numbers: an expected height for male recruits and an expected height for female recruits, each with a measure of uncertainty. That is a *prediction* question, posed as a quantity. + +### Exercise 4 + +In `analysis.qmd`, add a chunk that runs `predictions(fit_recruits)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r temperance-4} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r temperance-4-test} +#| echo: true +predictions(fit_recruits) +``` + +### + +`predictions()` runs the machine: it returns the expected outcome for each unit you give it, with uncertainty attached. With no `newdata`, it runs on every row of the data the model was fit on --- one prediction per recruit, so 50 rows for our 50-row sample. Each row's `estimate` is that recruit's model-predicted height. Notice already that the standard error column has just two distinct values --- about 1.01 for male rows and about 2.02 for female rows. We will come back to that. + +### + +One DGM answers many questions. Ask it for an expected value and it is a deterministic formula: same sex in, same number out --- which is why these 50 rows show only two distinct estimates. Ask it what will happen to a single recruit and the randomness matters: the DGM must draw from its residual distribution, and two recruits of the same sex can end up with very different heights. The question, not habit, picks the mode. + +### Exercise 5 + +Every marginaleffects call is either **marginal** --- run the machine on real recruits and average within groups --- or **conditional** --- run it on one representative, hypothetical recruit. The first uses `by`; the second uses `condition`. + +Add a chunk that runs `avg_predictions(fit_recruits, by = "sex")`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r temperance-5} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r temperance-5-test} +#| echo: true +avg_predictions(fit_recruits, by = "sex") +``` + +### + +`by` is always the same operation: run the model on every real recruit, then average within groups. This is a **marginal** prediction --- it uses the actual distribution of recruits in each group. The two rows it returns --- an average predicted height of about 166 cm for female recruits and about 175 cm for male recruits --- are the answer to the question we started with. + +### Exercise 6 + +Modify the chunk to run `plot_predictions(fit_recruits, condition = "sex")` instead. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r temperance-6} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r temperance-6-test} +#| echo: true +plot_predictions(fit_recruits, condition = "sex") +``` + +### + +With `condition`, the machine is run on imaginary recruits --- one representative recruit per sex, with any other variables held at typical values --- so the plot shows **conditional** expectations. Our model has only one covariate, so the conditional picture here nearly matches the marginal `avg_predictions(by = "sex")` numbers from the previous exercise; with more covariates the two can differ meaningfully. + +You should see two point estimates with 95% confidence intervals. The Female interval is visibly wider than the Male interval. That is the consequence of having drawn 10 female rows but 40 male rows: the standard error of a group mean shrinks roughly with the square root of the group's sample size, so quadrupling the sample size halves the standard error, and the female interval ends up about twice as wide. + +### Exercise 7 + +Modify the previous chunk to add `draw = FALSE` to the `plot_predictions()` call. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r temperance-7} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r temperance-7-test} +#| echo: true +plot_predictions(fit_recruits, condition = "sex", draw = FALSE) +``` + +### + +Because `plot_predictions()` returns a ggplot object, you can continue with ggplot commands like `labs()`. But it can be useful to see the underlying values in the tibble and build your own plot directly. + +### Exercise 8 + +Create a beautiful plot starting from the output of `plot_predictions(fit_recruits, condition = "sex", draw = FALSE)`. Do this in your QMD (much easier than typing in the R Terminal directly). Title: the key variables. Subtitle: the key takeaway. Caption: data source. Axis labels: human-readable with units. + +Paste the plot code below. + +```{r temperance-8} +question_text(NULL, + message = " plot_predictions(fit_recruits, condition = \"sex\", draw = FALSE) |>\n ggplot(aes(x = sex, y = estimate)) +\n geom_pointrange(aes(ymin = conf.low, ymax = conf.high),\n color = \"grey20\", size = 1.1) +\n labs(\n title = \"Forecast Height for USMC Recruits\",\n subtitle = \"Male recruits are expected to be about 9 cm taller\",\n x = \"Sex\",\n y = \"Expected Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 12) +``` + +### + +Here is our version: + +```{r} +#| echo: false +#| message: false +#| warning: false +plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |> + ggplot(aes(x = sex, y = estimate)) + + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), + color = "grey20", size = 1.1) + + labs( + title = "Forecast Height for USMC Recruits", + subtitle = "Male recruits are expected to be about 9 cm taller than female recruits", + x = "Sex", + y = "Expected Height (cm)", + caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" + ) + + theme_minimal() +``` + +Data science often involves this back-and-forth style of work. First, make a single chunk of code --- say, a new plot --- work well. This requires interactive work between the QMD and the R Terminal. Second, ensure that the entire QMD runs correctly on its own. + +### Exercise 9 + +Finalize the new graphics chunk in your QMD. Render to ensure it all works. In the R Terminal, run: + +``` +show_file("analysis.qmd", chunk = "Last") +``` + +CP/CR. + +```{r temperance-9} +question_text(NULL, + message = " > show_file(\"analysis.qmd\", chunk = \"Last\")\n ```{r}\n plot_predictions(fit_recruits, condition = \"sex\", draw = FALSE) |>\n ggplot(aes(x = sex, y = estimate)) +\n geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +\n labs(\n title = \"Forecast Height for USMC Recruits\",\n subtitle = \"Male recruits are expected to be about 9 cm taller\",\n x = \"Sex\", y = \"Expected Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()\n ```", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 8) +``` + +### + +The picture is only as trustworthy as the mechanism behind it. A beautiful graphic tells a story, but that story rests on assumptions --- validity, stability, representativeness --- that are never completely true. No `fit()` on the sample alone can rescue you from a misunderstood mechanism. + +### Exercise 10 + +Write the last sentence of your summary paragraph. It describes at least one Quantity of Interest and a measure of uncertainty. It is OK if this QoI differs from the one you began with. + +```{r temperance-10} +question_text(NULL, + message = "We estimate that the average male USMC recruit will be about 175 cm tall and the average female recruit about 166 cm, with the female interval visibly wider than the male interval because we have far fewer female observations in our sample.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Add a final sentence to your summary paragraph, but don't copy/paste our answer exactly. Render the QMD. + +### Exercise 11 + +Write a few sentences explaining why the estimates for the quantities of interest, and the uncertainty, might be wrong. Suggest alternative estimates and a wider confidence interval if warranted. + +```{r temperance-11} +question_text(NULL, + message = "USMC recruits are more physically fit and probably taller than the general young-adult population measured by NHANES. A better point estimate would nudge male recruit heights upward by a centimeter or two and female recruit heights similarly. The confidence intervals should be wider than what the model reports --- our 50-row sample is small (only 10 of those rows are female), and the bigger sources of uncertainty (the NHANES-to-USMC validity gap and the voluntary-participation representativeness gap) aren't captured by the model's standard error.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +Always go back to your Preceptor Table --- the information which, if you had it, would make answering your question easy. In almost all real-world cases, the Preceptor Table and the data are fairly different. So, even a perfectly estimated statistical model is rarely as useful as we might like. + +### Exercise 12 + +Rearrange the material in your QMD so the order is graphic, then summary paragraph. The chunk that creates `fit_recruits` must occur before the chunk that creates the graphic. You can keep or discard the math at your discretion. Render the QMD. In the R Terminal, run: + +``` +show_file("analysis.qmd") +``` + +CP/CR. + +```{r temperance-12} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 10) +``` + +### + +

+---
+title: "Recruits"
+author: "<your name>"
+execute:
+  echo: false
+---
+
+#| message: false
+library(tidyverse)
+library(tidymodels)
+library(primer.tutorials)
+library(broom)
+library(marginaleffects)
+
+#| cache: true
+fit_recruits <- linear_reg() |>
+  set_engine("lm") |>
+  fit(height ~ sex, data = recruits)
+
+plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |>
+  ggplot(aes(x = sex, y = estimate)) +
+  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
+  labs(
+    title = "Forecast Height for USMC Recruits",
+    subtitle = "Male recruits are expected to be about 9 cm taller",
+    x = "Sex", y = "Expected Height (cm)",
+    caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits"
+  ) +
+  theme_minimal()
+
+People vary in height, and that variation is strongly patterned by sex. Using a 50-row teaching sample drawn from the National Health and Nutrition Examination Survey (NHANES, conducted by the Centers for Disease Control and Prevention), we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. We model height as a normally distributed variable which is a linear function of sex. We estimate that the average male recruit will be about 175 cm tall and the average female recruit about 166 cm, with confidence intervals that are visibly wider for the female estimate because there are only ten female rows in our sample.
+
+ +### + +This is the version of your QMD file at which your teacher is most likely to look closely. + +### Exercise 13 + +Publish your rendered QMD to GitHub Pages. In the bash Terminal (not the R Terminal!), run: + +``` +quarto publish gh-pages analysis.qmd +``` + +Copy/paste the resulting URL below. + +```{r temperance-13} +question_text(NULL, + message = "https://.github.io/recruits/analysis.html", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 2) +``` + +### + +Commit/push everything. + +### Exercise 14 + +Copy/paste the URL to your GitHub repo. + +```{r temperance-14} +question_text(NULL, + message = "https://github.com//recruits", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 2) +``` + +### + +A fitted DGM is a machine for generating answers, and using it well means asking well: pose the question as a quantity, specify the units --- real or imaginary --- whose outcomes you care about, then run the machine and draw the picture. If all our assumptions were correct, our DGM would be true, and there would be no better way to predict the future. Sadly, that is only ever the case with toy examples involving coins and dice. + +The world confronts us. Make decisions we must. + + +## Summary +### + +```{r} +#| echo: false +#| message: false +#| warning: false +plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |> + ggplot(aes(x = sex, y = estimate)) + + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), + color = "grey20", size = 1.1) + + labs( + title = "Forecast Height for USMC Recruits", + subtitle = "Male recruits are expected to be about 9 cm taller than female recruits", + x = "Sex", + y = "Expected Height (cm)", + caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" + ) + + theme_minimal() +``` + +People vary in height, and that variation is strongly patterned by sex. Using a 50-row teaching sample drawn from the National Health and Nutrition Examination Survey (NHANES, conducted by the Centers for Disease Control and Prevention), we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. We model height as a normally distributed variable which is a linear function of sex. We estimate that the average male recruit will be about 175 cm tall and the average female recruit about 166 cm, with the female interval visibly wider than the male interval because we have only ten female rows in our sample. + +This forecast would help you estimate how many small-versus-large uniforms to order. It does not address every question the logistics office will face --- it says nothing about weight, shoulder width, or inseam. But it gives a defensible starting point. + +The world is always more uncertain than our models would have us believe. + +```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} +``` diff --git a/inst/tutorials/05-recruits-wisdom/tutorial.Rmd b/inst/tutorials/05-recruits-wisdom/tutorial.Rmd new file mode 100644 index 0000000..b868ee8 --- /dev/null +++ b/inst/tutorials/05-recruits-wisdom/tutorial.Rmd @@ -0,0 +1,449 @@ +--- +title: "Recruits: Wisdom" +author: David Kane +tutorial: + id: 05-recruits-wisdom +output: + learnr::tutorial: + progressive: yes + allow_skip: yes +runtime: shiny_prerendered +description: "Tutorial #05 (Wisdom) for Preceptor's Primer" +--- + +```{r setup, include = FALSE} +library(learnr) +library(tutorial.helpers) +library(gt) + +library(tidyverse) +library(primer.tutorials) + +knitr::opts_chunk$set(echo = FALSE) +knitr::opts_chunk$set(out.width = "90%") +options(tutorial.exercise.timelimit = 600, + tutorial.storage = "local") +``` + +```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} +``` + + +## Introduction +### + +This tutorial supports the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science: Using the Cardinal Virtues for Inference*](https://ppbds.github.io/primer/) by [David Kane](https://davidkane.info/). + +The world confronts us. Make decisions we must. + +Imagine that you are a logistics analyst at U.S. Marine Corps Supply, planning next year's bootcamp uniform order. The Quartermaster General has a big-picture goal --- every recruit kitted out at the right size on day one, within a fixed procurement budget --- and trusts you to figure out what to order and in what quantities. Many models would inform the order: distributions of height, weight, chest, foot length; the male-female mix in incoming recruits; the size-tail dynamics that drive how many extra-large uniforms to keep in inventory. This tutorial builds just one of them: the average height of male and female recruits. The number alone won't determine the order, but the height-by-sex split is one good input to the larger forecast. There are many decisions to make. + +### + +This is the first of four Recruits tutorials, one per Cardinal Virtue: **Wisdom** (the question, the Preceptor Table, and a first look at the data), **Justice** (an audit of assumptions and the choice of model family), **Courage** (building the data generating mechanism), and **Temperance** (the answer, with its uncertainty). You will build one `analysis.qmd` across all four. + +### Exercise 1 + +You should be connected to a repo named `recruits`. If you are not, create one and connect to it now. Create a new Quarto document titled `"Recruits"` with yourself as the author, save it as `analysis.qmd`, render it, and open `analysis.html` with Live Server so the rendered HTML auto-refreshes on every later render. + +Put `library(tidyverse)` and `library(primer.tutorials)` in a setup chunk with `#| message: false`, and add + +``` +execute: + echo: false +``` + +to the YAML header. Create a `.gitignore` file with `analysis_files` on the first line followed by a blank line. Commit and push. Use AI however you like. + +In the R Terminal, run `show_file(".gitignore")`. If that fails, it is probably because you have not yet loaded `library(tutorial.helpers)` in the R Terminal. + +"CP/CR" is short for *copy/paste the **c**ommand and the **r**esponse*. CP/CR. + +```{r introduction-1} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +Rendering a QMD converts the code and prose into an HTML document that shows only what you want readers to see --- in this case, the eventual plot and summary paragraph, with the code hidden. Professional data-science deliverables almost never show raw code; the code is the path, not the destination. + +### Exercise 2 + +Add a code chunk to `analysis.qmd` that prints the `recruits` tibble (just `recruits` on its own line). Render, and copy the first few rows from the rendered HTML below. + +```{r introduction-2} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 6) +``` + +### + +The `recruits` tibble is a 50-row teaching cut of NHANES --- the National Health and Nutrition Examination Survey, conducted by the Centers for Disease Control and Prevention --- restricted to ages 18--27 and to the columns `height` (centimeters), `sex`, and `age`. It has 40 male and 10 female rows. That split is deliberately uneven, not a reflection of the real male/female mix among recruits; it will matter for the precision of the two group means in Temperance. + +NHANES also intentionally oversamples some demographic groups --- older adults, racial and ethnic minorities --- so each subgroup has enough observations to study on its own. The raw survey is not representative of the general population; NHANES uses weighting to correct for this. Our `recruits` sample inherits and amplifies these issues, which we flag in Justice. + + +## Wisdom +### + +Data science starts with some broad questions and a data set which might help us to answer them. + +The specific question: *What is the average height of male and female USMC recruits?* + +We will work from the `recruits` tibble in the **[primer.tutorials](https://github.com/PPBDS/primer.tutorials)** package --- a teaching cut of the National Health and Nutrition Examination Survey (NHANES), conducted by the Centers for Disease Control and Prevention, drawn from the `nhanes` tibble of the **[primer.data](https://github.com/PPBDS/primer.data)** package. + +### Exercise 1 + +In your own words, describe the key components of Wisdom when working on a data science problem. + +```{r wisdom-1} +question_text(NULL, + message = "Wisdom begins with a question and then moves on to the creation of a Preceptor Table and an examination of our data.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +> *The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.* --- John W. Tukey + +### Exercise 2 + +Define a Preceptor Table. + +```{r wisdom-2} +question_text(NULL, + message = "A Preceptor Table is the smallest possible table of data with rows and columns such that, if there is no missing data, we can easily calculate the quantity of interest.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +The Preceptor Table does not include all the covariates which you will eventually include in your model. It only includes, along with the outcome(s), covariates which are mentioned in your question. + +### Exercise 3 + +Our first look at the data always starts with the outcome variable. Write R code that plots the distribution of `height` in `recruits` as a density plot, with a descriptive title, subtitle, axis labels, and caption. Paste the code below and run it. + +```{r wisdom-3} +question_text(NULL, + message = " recruits |>\n ggplot(aes(x = height)) +\n geom_density(fill = \"grey70\", color = \"grey30\") +\n labs(\n title = \"Height in Our 50-Recruit Sample\",\n subtitle = \"A broad, asymmetric peak with a heavy left shoulder --- a hint that sex matters for height\",\n x = \"Height (cm)\", y = NULL,\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal() +\n theme(axis.text.y = element_blank(),\n axis.ticks.y = element_blank())", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 10) +``` + +### + +```{r} +#| echo: false +#| message: false +#| warning: false +recruits |> + ggplot(aes(x = height)) + + geom_density(fill = "grey70", color = "grey30") + + labs( + title = "Height in Our 50-Recruit Sample", + subtitle = "A broad, asymmetric peak with a heavy left shoulder --- a hint that sex matters for height", + x = "Height (cm)", + y = NULL, + caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" + ) + + theme_minimal() + + theme(axis.text.y = element_blank(), + axis.ticks.y = element_blank()) +``` + +### Exercise 4 + +Now look at how `height` varies with `sex`. Write R code that plots the individual `height` values in `recruits` as a jitter plot, with `sex` on the x-axis and `height` on the y-axis, with a descriptive title, subtitle, axis labels, and caption. Paste the code below and run it. + +```{r wisdom-4} +question_text(NULL, + message = " recruits |>\n ggplot(aes(x = sex, y = height)) +\n geom_jitter(width = 0.15, alpha = 0.5) +\n labs(\n title = \"Height by Sex in Our 50-Recruit Sample\",\n subtitle = \"Male recruits average about 9 cm taller, and there are far fewer women\",\n x = \"Sex\", y = \"Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 10) +``` + +### + +```{r} +#| echo: false +#| message: false +#| warning: false +recruits |> + ggplot(aes(x = sex, y = height)) + + geom_jitter(width = 0.15, alpha = 0.5) + + labs( + title = "Height by Sex in Our 50-Recruit Sample", + subtitle = "Male recruits average about 9 cm taller, and there are far fewer women", + x = "Sex", + y = "Height (cm)", + caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" + ) + + theme_minimal() +``` + +### Exercise 5 + +Describe the key components of Preceptor Tables in general, without worrying about this specific problem. Use words like "units," "outcomes," and "covariates." + +```{r wisdom-5} +question_text(NULL, + message = "The rows of the Preceptor Table are the units. The outcome is at least one of the columns. If the problem is causal, there will be at least two (potential) outcome columns. The other columns are covariates. If the problem is causal, at least one of the covariates will be considered a treatment.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 6) +``` + +### + +This problem is predictive, so there is only one outcome column: `height`. We compare average heights between male and female young adults. + +### Exercise 6 + +What are the units for this problem? + +```{r wisdom-6} +question_text(NULL, + message = "Each row represents one USMC recruit expected to enlist during the coming year. There are roughly 5,000 such recruits per year.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +Specifying the Preceptor Table forces us to think clearly about the units and outcomes implied by the question. No data science project follows a single direction. We always backtrack. There is always dialogue. We model units, but we only really care about aggregates. + +### Exercise 7 + +What is the outcome variable for this problem? + +```{r wisdom-7} +question_text(NULL, + message = "Keep track of two 'outcome' variables: the one in our Preceptor Table and the one in our data. In this case, they share the name `height` --- both refer to a young adult's standing height in centimeters --- though the *measurement procedure* behind each column differs.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +The outcome variable that we really care about is often not the outcome variable which our data includes. This compromise --- working with what we have rather than what we really want --- is a part of most data science work in the real world. + +### Exercise 8 + +What is a covariate which you think might be useful for this problem, regardless of whether or not it might be included in the data? + +```{r wisdom-8} +question_text(NULL, + message = "Age at enlistment is a candidate. Height is nearly-but-not-quite stable by the age-18 window we care about, so a recruit ten years older is almost the same height as a fresh one. Another candidate is race or ethnicity, since height distributions vary across ancestry groups.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +The term "covariates" is used in at least three ways in data science. First, it is all the variables which *might* be useful, regardless of whether or not we have the data. Second, it is all the variables for which we have data. Third, it is the set of variables in the data which we end up using in the model. + +### Exercise 9 + +What are the treatments, if any, for this problem? + +```{r wisdom-9} +question_text(NULL, + message = "There are no treatments in this problem because we are building a predictive model. Nothing in `recruits` can be manipulated --- we cannot assign a recruit a different sex or a different height. Predictive models have covariates, not treatments.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +Remember that a treatment is just another covariate which, for the purposes of this specific problem, we are assuming can be manipulated --- thereby creating two or more different potential outcomes for each unit. Predictive models have no such variable. + +### Exercise 10 + +What moment in time does the Preceptor Table refer to? + +```{r wisdom-10} +question_text(NULL, + message = "The coming year --- specifically, the moment at which each recruit arrives at bootcamp and has their height recorded at enlistment.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 3) +``` + +### + +A Preceptor Table can never really refer to an exact instant in time since nothing is instantaneous in this fallen world. In practice, the Preceptor Table's "moment" is a narrow window: recruit intake for a single calendar year. + +### Exercise 11 + +Describe in words the Preceptor Table for this problem. + +```{r wisdom-11} +question_text(NULL, + message = "The Preceptor Table has three columns: one identifying each USMC recruit, one for `height`, and one for `sex`. Each row represents one recruit expected to enlist during the coming year.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +```{r} +#| echo: false +#| results: asis +pre_tibble <- tibble::tribble( + ~`Recruit`, ~`Height (cm)`, ~`Sex`, + "Maya Castillo", "150", "Female", + "Tyler Brooks", "170", "Male", + "...", "...", "...", + "Grace Sullivan", "160", "Female" +) + +pre_title_footnote <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." +pre_units_footnote <- "Each row is one USMC recruit expected to enlist during the coming year. The example rows use plausible names; missing rows represent the other roughly 5,000 recruits not shown." +pre_outcome_footnote <- "Height in centimeters, measured at enlistment." +pre_covariate_footnote <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." + +pre_gt_html <- gt::gt(pre_tibble, id = "preceptor_tbl") |> + gt::tab_header(title = "Preceptor Table") |> + gt::tab_spanner(label = "Unit", id = "unit_span", + columns = c(`Recruit`)) |> + gt::tab_spanner(label = "Outcome", id = "outcome_span", + columns = c(`Height (cm)`)) |> + gt::tab_spanner(label = "Covariate", id = "covariates_span", + columns = c(`Sex`)) |> + gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> + gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> + gt::fmt_markdown(columns = gt::everything()) |> + gt::tab_options( + heading.title.font.size = "1.5em", + heading.title.font.weight = "bolder", + column_labels.font.weight = "normal" + ) |> + gt::tab_style( + style = gt::cell_text(weight = "bold"), + locations = gt::cells_column_spanners() + ) |> + gt::tab_style( + style = gt::cell_text(align = "left"), + locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) + ) |> + gt::tab_style( + style = gt::cell_text(align = "right"), + locations = gt::cells_column_spanners(spanners = "outcome_span") + ) |> + gt::tab_footnote(footnote = pre_title_footnote, + locations = gt::cells_title()) |> + gt::tab_footnote(footnote = pre_units_footnote, + locations = gt::cells_column_spanners(spanners = "unit_span")) |> + gt::tab_footnote(footnote = pre_outcome_footnote, + locations = gt::cells_column_spanners(spanners = "outcome_span")) |> + gt::tab_footnote(footnote = pre_covariate_footnote, + locations = gt::cells_column_spanners(spanners = "covariates_span")) |> + gt::opt_css( + css = " + #preceptor_tbl .gt_footnote { + max-width: 1px; + word-break: break-word; + } + " + ) |> + gt::as_raw_html(inline_css = FALSE) + +cat( + "```{=html}\n", + '
', + pre_gt_html, + "
\n", + "```\n", + sep = "" +) +``` + +### Exercise 12 + +Because `recruits` is already prepared (filtered to ages 18--27, restricted to `height`, `sex`, and `age`, and missing rows dropped), there is no data-prep code to write. Take a moment to verify the breakdown that matters most for what comes next. In `analysis.qmd`, add a chunk that runs `recruits |> count(sex)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. + +```{r wisdom-12} +question_text(NULL, + answer(NULL, correct = TRUE), + allow_retry = TRUE, + try_again_button = "Edit Answer", + incorrect = NULL, + rows = 5) +``` + +### + +```{r wisdom-12-test} +#| echo: true +recruits |> count(sex) +``` + +### + +The split between sexes is deliberate: 40 male rows and 10 female. We did not choose this to reflect the male/female proportions among real USMC recruits, but to make a particular feature of the model visible later --- when we estimate the average height for each group, the larger group will give a tighter estimate (smaller confidence interval) than the smaller group. + +> *You can never look at the data too much.* --- Mark Engerman + +### Exercise 13 + +We will be creating a summary paragraph over the course of these tutorials. Write the first two sentences. The first sentence is a general statement about the overall topic, mentioning the general class of outcome and at least one covariate. The second sentence introduces the data source and the specific question --- when/where gathered, how many observations, who collected it. + +```{r wisdom-13} +question_text(NULL, + message = "People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year.", + answer(NULL, correct = TRUE), + allow_retry = FALSE, + incorrect = NULL, + rows = 5) +``` + +### + +Read our answer. It will not be the same as yours. You can change your answer to incorporate some of our ideas, but do not copy/paste our answer exactly. Add your two sentences to the summary paragraph portion of your QMD. Render the QMD, then commit and push. + + +## Summary +### + +This tutorial covered the **Wisdom** stage of the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science*](https://ppbds.github.io/primer/). + +Wisdom begins with a question --- *What is the average height of male and female USMC recruits?* --- and then builds a Preceptor Table and examines the data. The Preceptor Table for this problem has one row per recruit expected to enlist in the coming year and three columns: the recruit, `height`, and `sex`. A first look at the `recruits` data --- a density plot of `height` and a jitter plot of `height` by `sex` --- shows male recruits averaging about 9 cm taller, with far fewer women in the sample. + +The **Justice** tutorial audits the assumptions behind this setup and chooses the model's functional form. + +The world is always more uncertain than our models would have us believe. + +```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} +``` diff --git a/inst/tutorials/05-recruits/tutorial.Rmd b/inst/tutorials/05-recruits/tutorial.Rmd deleted file mode 100644 index d9b2b93..0000000 --- a/inst/tutorials/05-recruits/tutorial.Rmd +++ /dev/null @@ -1,1812 +0,0 @@ ---- -title: Recruits -author: David Kane -tutorial: - id: 05-recruits -output: - learnr::tutorial: - progressive: yes - allow_skip: yes -runtime: shiny_prerendered -description: "Tutorial #05 for Preceptor's Primer" ---- - -```{r setup, include = FALSE} -library(learnr) -library(tutorial.helpers) -library(gt) - -library(tidyverse) -library(tidymodels) -library(primer.tutorials) -library(broom) -library(marginaleffects) - -knitr::opts_chunk$set(echo = FALSE) -knitr::opts_chunk$set(out.width = "90%") -options(tutorial.exercise.timelimit = 600, - tutorial.storage = "local") - -fit_recruits <- linear_reg() |> - set_engine("lm") |> - fit(height ~ sex, data = recruits) -``` - -```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} -``` - - -## Introduction -### - -This tutorial supports the [Recruits](https://ppbds.github.io/primer/05-recruits.html) chapter of [*Preceptor's Primer for Bayesian Data Science: Using the Cardinal Virtues for Inference*](https://ppbds.github.io/primer/) by [David Kane](https://davidkane.info/). - -The world confronts us. Make decisions we must. - -Imagine that you are a logistics analyst at U.S. Marine Corps Supply, planning next year's bootcamp uniform order. The Quartermaster General has a big-picture goal --- every recruit kitted out at the right size on day one, within a fixed procurement budget --- and trusts you to figure out what to order and in what quantities. Many models would inform the order: distributions of height, weight, chest, foot length; the male-female mix in incoming recruits; the size-tail dynamics that drive how many extra-large uniforms to keep in inventory. This tutorial builds just one of them: the average height of male and female recruits. The number alone won't determine the order, but the height-by-sex split is one good input to the larger forecast. There are many decisions to make. - -### Exercise 1 - -What are the four [Cardinal Virtues](https://en.wikipedia.org/wiki/Cardinal_virtues), in order, which we use to guide our data science work? - -```{r introduction-1} -question_text(NULL, - message = "Wisdom, Justice, Courage, and Temperance.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 2) -``` - -### - -The four virtues are the four sections of every Primer chapter and tutorial: a question and a Preceptor Table (Wisdom), an audit of assumptions (Justice), the model (Courage), and the answer with its uncertainty (Temperance). *Recruits* is the first time we run that arc end to end, on a 50-row sample of height by sex. - -### Exercise 2 - -You should be working inside a GitHub repo named `recruits`, opened in a Codespace from the [`PPBDS/codespace-starter`](https://github.com/PPBDS/codespace-starter) template. If you are not there yet, please create that repo and open it in a Codespace now --- see the [package README](https://github.com/PPBDS/primer/tree/main/primer.tutorials#working-environments-and-repo-setup) if you need setup instructions or are working locally instead. - -Once you're inside the `recruits` repo, create a new Quarto document titled `"Recruits"` with yourself as the author, render it, and save it as `analysis.qmd`. - -Create a `.gitignore` file with `analysis_files` on the first line followed by a blank line. Commit and push. - -In the R Terminal, run: - -``` -show_file(".gitignore") -``` - -If that fails, it's probably because you haven't loaded `library(tutorial.helpers)` in the R Terminal. - -"CP/CR" is short for *copy/paste the **c**ommand and the **r**esponse* --- both the command you sent to R and the response R gave back. - -CP/CR. - -```{r introduction-2} -question_text(NULL, - message = " > show_file(\".gitignore\")\n analysis_files", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -The `recruits` tibble is a 50-row teaching cut of NHANES with 40 male and 10 female young adults aged 18--27; the deliberately uneven split will matter when we get to the precision of group means in Temperance. - -### Exercise 3 - -In your QMD, put `library(tidyverse)` and `library(primer.tutorials)` in a new code chunk. Render the file. - -Notice that the file does not look good because the code is visible and there are messages cluttering the output. To fix this, add `#| message: false` to this setup chunk to remove the messages. Also add the following to the YAML header, between the `---` delimiters, to hide the code itself: - -``` -execute: - echo: false -``` - -Render again --- it should look cleaner. In the R Terminal, run: - -``` -show_file("analysis.qmd", chunk = "Last") -``` - -CP/CR. - -```{r introduction-3} -question_text(NULL, - message = " > show_file(\"analysis.qmd\", chunk = \"Last\")\n ```{r}\n #| message: false\n library(tidyverse)\n library(primer.tutorials)\n ```", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Rendering a QMD converts the code and prose into an HTML document that shows only what you want readers to see --- in this case, the eventual plot and summary paragraph, with the code hidden. Professional data-science deliverables almost never show raw code; the code is the path, not the destination. - -### Exercise 4 - -In the R Terminal, type `?recruits` and paste the Description below. - -```{r introduction-4} -question_text(NULL, - message = "A 50-row sample drawn from the National Health and Nutrition Examination Survey (NHANES, via primer.data::nhanes) for the 05-recruits tutorial. Restricted to ages 18--27 and to the columns relevant for the tutorial's question. The split between sexes (40 male, 10 female) is deliberately uneven so the two group means have visibly different standard errors --- a feature the tutorial's Temperance section asks students to notice and explain.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Note that NHANES intentionally oversamples some demographic groups --- older adults, racial and ethnic minorities --- so that each subgroup has enough observations to be studied on its own. The raw NHANES survey is not representative of the general population; the survey uses weighting to correct for this. - -Our `recruits` sample inherits these issues and amplifies them: we drew 40 men and 10 women rather than matching their representation in NHANES. We should flag both as concerns later in Justice. - -### Exercise 5 - -Height of young adults is the broad topic of this tutorial. Given that topic, which variable in `recruits` should we use as our outcome variable? - -```{r introduction-5} -question_text(NULL, - message = "Our outcome variable is `height`. We care about the height of USMC recruits, who are young adults.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -We will use `height` as our outcome variable. - -### Exercise 6 - -Which variable in `recruits` do you think might have an important connection to `height`? - -```{r introduction-6} -question_text(NULL, - message = "One candidate is `sex`. Men and women tend to have systematically different average heights, so a model that includes `sex` should do better than one that does not.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -With a predictive model, each individual unit has only one observed outcome. Predictive models have no "treatments" --- only covariates. - -### Exercise 7 - -Specify two different groups of young adults which have different values for `sex` and which might have different average values for `height`. - -```{r introduction-7} -question_text(NULL, - message = "Consider two groups: young adults with `sex` = Male, and young adults with `sex` = Female. These two groups might have different average heights.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -In predictive models, do not use "cause," "influence," "impact," or anything else which suggests causation. The best phrasing is in terms of "differences" between groups of units with different values for a covariate of interest. - -### Exercise 8 - -Write a predictive question connecting the outcome `height` to `sex`, the covariate of interest. - -```{r introduction-8} -question_text(NULL, - message = "What is the average height of male and female USMC recruits?", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -The average height by sex is a convenient question --- a single number per sex, easy to interpret, easy to use for ordering uniforms. It isn't the only useful question. A logistics officer might also want the tallest recruit expected next year (so no one arrives without a uniform that fits), or the 10th and 90th percentiles of each sex's height distribution (to know how many small-versus-large sizes to order). We stick with the average question here because it exercises every step of the Cardinal Virtues cleanly; the matching chapter considers the broader family. - - -## Wisdom -### - -Data science starts with some broad questions and a data set which might help us to answer them. - -The specific question: *What is the average height of male and female USMC recruits?* - -We will work from the `recruits` tibble in the **[primer.tutorials](https://github.com/PPBDS/primer.tutorials)** package --- a teaching cut of the National Health and Nutrition Examination Survey (NHANES), conducted by the Centers for Disease Control and Prevention, drawn from the `nhanes` tibble of the **[primer.data](https://github.com/PPBDS/primer.data)** package. - -### Exercise 1 - -In your own words, describe the key components of Wisdom when working on a data science problem. - -```{r wisdom-1} -question_text(NULL, - message = "Wisdom begins with a question and then moves on to the creation of a Preceptor Table and an examination of our data.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -> *The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.* --- John W. Tukey - -### Exercise 2 - -Define a Preceptor Table. - -```{r wisdom-2} -question_text(NULL, - message = "A Preceptor Table is the smallest possible table of data with rows and columns such that, if there is no missing data, we can easily calculate the quantity of interest.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -The Preceptor Table does not include all the covariates which you will eventually include in your model. It only includes, along with the outcome(s), covariates which are mentioned in your question. - -### Exercise 3 - -Our first look at the data always starts with the outcome variable. Write R code that plots the distribution of `height` in `recruits` as a density plot, with a descriptive title, subtitle, axis labels, and caption. Paste the code below and run it. - -```{r wisdom-3} -question_text(NULL, - message = " recruits |>\n ggplot(aes(x = height)) +\n geom_density(fill = \"grey70\", color = \"grey30\") +\n labs(\n title = \"Height in Our 50-Recruit Sample\",\n subtitle = \"A broad, asymmetric peak with a heavy left shoulder --- a hint that sex matters for height\",\n x = \"Height (cm)\", y = NULL,\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal() +\n theme(axis.text.y = element_blank(),\n axis.ticks.y = element_blank())", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 10) -``` - -### - -```{r} -#| echo: false -#| message: false -#| warning: false -recruits |> - ggplot(aes(x = height)) + - geom_density(fill = "grey70", color = "grey30") + - labs( - title = "Height in Our 50-Recruit Sample", - subtitle = "A broad, asymmetric peak with a heavy left shoulder --- a hint that sex matters for height", - x = "Height (cm)", - y = NULL, - caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" - ) + - theme_minimal() + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank()) -``` - -### Exercise 4 - -Now look at how `height` varies with `sex`. Write R code that plots the individual `height` values in `recruits` as a jitter plot, with `sex` on the x-axis and `height` on the y-axis, with a descriptive title, subtitle, axis labels, and caption. Paste the code below and run it. - -```{r wisdom-4} -question_text(NULL, - message = " recruits |>\n ggplot(aes(x = sex, y = height)) +\n geom_jitter(width = 0.15, alpha = 0.5) +\n labs(\n title = \"Height by Sex in Our 50-Recruit Sample\",\n subtitle = \"Male recruits average about 9 cm taller, and there are far fewer women\",\n x = \"Sex\", y = \"Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 10) -``` - -### - -```{r} -#| echo: false -#| message: false -#| warning: false -recruits |> - ggplot(aes(x = sex, y = height)) + - geom_jitter(width = 0.15, alpha = 0.5) + - labs( - title = "Height by Sex in Our 50-Recruit Sample", - subtitle = "Male recruits average about 9 cm taller, and there are far fewer women", - x = "Sex", - y = "Height (cm)", - caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" - ) + - theme_minimal() -``` - -### Exercise 5 - -Describe the key components of Preceptor Tables in general, without worrying about this specific problem. Use words like "units," "outcomes," and "covariates." - -```{r wisdom-5} -question_text(NULL, - message = "The rows of the Preceptor Table are the units. The outcome is at least one of the columns. If the problem is causal, there will be at least two (potential) outcome columns. The other columns are covariates. If the problem is causal, at least one of the covariates will be considered a treatment.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -This problem is predictive, so there is only one outcome column: `height`. We compare average heights between male and female young adults. - -### Exercise 6 - -What are the units for this problem? - -```{r wisdom-6} -question_text(NULL, - message = "Each row represents one USMC recruit expected to enlist during the coming year. There are roughly 5,000 such recruits per year.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -Specifying the Preceptor Table forces us to think clearly about the units and outcomes implied by the question. No data science project follows a single direction. We always backtrack. There is always dialogue. We model units, but we only really care about aggregates. - -### Exercise 7 - -What is the outcome variable for this problem? - -```{r wisdom-7} -question_text(NULL, - message = "Keep track of two 'outcome' variables: the one in our Preceptor Table and the one in our data. In this case, they share the name `height` --- both refer to a young adult's standing height in centimeters --- though the *measurement procedure* behind each column differs.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -The outcome variable that we really care about is often not the outcome variable which our data includes. This compromise --- working with what we have rather than what we really want --- is a part of most data science work in the real world. - -### Exercise 8 - -What is a covariate which you think might be useful for this problem, regardless of whether or not it might be included in the data? - -```{r wisdom-8} -question_text(NULL, - message = "Age at enlistment is a candidate. Height is nearly-but-not-quite stable by the age-18 window we care about, so a recruit ten years older is almost the same height as a fresh one. Another candidate is race or ethnicity, since height distributions vary across ancestry groups.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -The term "covariates" is used in at least three ways in data science. First, it is all the variables which *might* be useful, regardless of whether or not we have the data. Second, it is all the variables for which we have data. Third, it is the set of variables in the data which we end up using in the model. - -### Exercise 9 - -What are the treatments, if any, for this problem? - -```{r wisdom-9} -question_text(NULL, - message = "There are no treatments in this problem because we are building a predictive model. Nothing in `recruits` can be manipulated --- we cannot assign a recruit a different sex or a different height. Predictive models have covariates, not treatments.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -Remember that a treatment is just another covariate which, for the purposes of this specific problem, we are assuming can be manipulated --- thereby creating two or more different potential outcomes for each unit. Predictive models have no such variable. - -### Exercise 10 - -What moment in time does the Preceptor Table refer to? - -```{r wisdom-10} -question_text(NULL, - message = "The coming year --- specifically, the moment at which each recruit arrives at bootcamp and has their height recorded at enlistment.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -A Preceptor Table can never really refer to an exact instant in time since nothing is instantaneous in this fallen world. In practice, the Preceptor Table's "moment" is a narrow window: recruit intake for a single calendar year. - -### Exercise 11 - -Describe in words the Preceptor Table for this problem. - -```{r wisdom-11} -question_text(NULL, - message = "The Preceptor Table has three columns: one identifying each USMC recruit, one for `height`, and one for `sex`. Each row represents one recruit expected to enlist during the coming year.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -```{r} -#| echo: false -#| results: asis -pre_tibble <- tibble::tribble( - ~`Recruit`, ~`Height (cm)`, ~`Sex`, - "Maya Castillo", "150", "Female", - "Tyler Brooks", "170", "Male", - "...", "...", "...", - "Grace Sullivan", "160", "Female" -) - -pre_title_footnote <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." -pre_units_footnote <- "Each row is one USMC recruit expected to enlist during the coming year. The example rows use plausible names; missing rows represent the other roughly 5,000 recruits not shown." -pre_outcome_footnote <- "Height in centimeters, measured at enlistment." -pre_covariate_footnote <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." - -pre_gt_html <- gt::gt(pre_tibble, id = "preceptor_tbl") |> - gt::tab_header(title = "Preceptor Table") |> - gt::tab_spanner(label = "Unit", id = "unit_span", - columns = c(`Recruit`)) |> - gt::tab_spanner(label = "Outcome", id = "outcome_span", - columns = c(`Height (cm)`)) |> - gt::tab_spanner(label = "Covariate", id = "covariates_span", - columns = c(`Sex`)) |> - gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> - gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> - gt::fmt_markdown(columns = gt::everything()) |> - gt::tab_options( - heading.title.font.size = "1.5em", - heading.title.font.weight = "bolder", - column_labels.font.weight = "normal" - ) |> - gt::tab_style( - style = gt::cell_text(weight = "bold"), - locations = gt::cells_column_spanners() - ) |> - gt::tab_style( - style = gt::cell_text(align = "left"), - locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) - ) |> - gt::tab_style( - style = gt::cell_text(align = "right"), - locations = gt::cells_column_spanners(spanners = "outcome_span") - ) |> - gt::tab_footnote(footnote = pre_title_footnote, - locations = gt::cells_title()) |> - gt::tab_footnote(footnote = pre_units_footnote, - locations = gt::cells_column_spanners(spanners = "unit_span")) |> - gt::tab_footnote(footnote = pre_outcome_footnote, - locations = gt::cells_column_spanners(spanners = "outcome_span")) |> - gt::tab_footnote(footnote = pre_covariate_footnote, - locations = gt::cells_column_spanners(spanners = "covariates_span")) |> - gt::opt_css( - css = " - #preceptor_tbl .gt_footnote { - max-width: 1px; - word-break: break-word; - } - " - ) |> - gt::as_raw_html(inline_css = FALSE) - -cat( - "```{=html}\n", - '
', - pre_gt_html, - "
\n", - "```\n", - sep = "" -) -``` - -### Exercise 12 - -Because `recruits` is already prepared (filtered to ages 18--27, restricted to `height`, `sex`, and `age`, and missing rows dropped), there is no data-prep code to write. Take a moment to verify the breakdown that matters most for what comes next. In `analysis.qmd`, add a chunk that runs `recruits |> count(sex)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r wisdom-12} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r wisdom-12-test} -#| echo: true -recruits |> count(sex) -``` - -### - -The split between sexes is deliberate: 40 male rows and 10 female. We did not choose this to reflect the male/female proportions among real USMC recruits, but to make a particular feature of the model visible later --- when we estimate the average height for each group, the larger group will give a tighter estimate (smaller confidence interval) than the smaller group. - -> *You can never look at the data too much.* --- Mark Engerman - -### Exercise 13 - -We will be creating a summary paragraph over the course of this tutorial. Write the first two sentences. The first sentence is a general statement about the overall topic, mentioning the general class of outcome and at least one covariate. The second sentence introduces the data source and the specific question --- when/where gathered, how many observations, who collected it. - -```{r wisdom-13} -question_text(NULL, - message = "People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Read our answer. It will not be the same as yours. You can change your answer to incorporate some of our ideas, but do not copy/paste our answer exactly. Add your two sentences to the summary paragraph portion of your QMD. Render the QMD. - - - -## Justice -### - -Wisdom gives us the Preceptor Table and the data. - -```{r} -#| echo: false -#| results: asis -pre_tibble_j <- tibble::tribble( - ~`Recruit`, ~`Height (cm)`, ~`Sex`, - "Maya Castillo", "150", "Female", - "Tyler Brooks", "170", "Male", - "...", "...", "...", - "Grace Sullivan", "160", "Female" -) - -pre_title_footnote_j <- "A Preceptor Table is the smallest table such that, if every cell were filled with its true value, answering the question 'What is the average height of male and female USMC recruits?' would be easy." -pre_units_footnote_j <- "Each row is one USMC recruit expected to enlist during the coming year. Missing rows represent the other recruits not shown." -pre_outcome_footnote_j <- "Height in centimeters, measured at enlistment." -pre_covariate_footnote_j <- "Sex, taking values 'Male' or 'Female', as recorded at enlistment." - -pre_gt_html_j <- gt::gt(pre_tibble_j, id = "preceptor_tbl_justice") |> - gt::tab_header(title = "Preceptor Table") |> - gt::tab_spanner(label = "Unit", id = "unit_span", - columns = c(`Recruit`)) |> - gt::tab_spanner(label = "Outcome", id = "outcome_span", - columns = c(`Height (cm)`)) |> - gt::tab_spanner(label = "Covariate", id = "covariates_span", - columns = c(`Sex`)) |> - gt::cols_align(align = "left" , columns = c(`Recruit`, `Sex`)) |> - gt::cols_align(align = "right", columns = c(`Height (cm)`)) |> - gt::fmt_markdown(columns = gt::everything()) |> - gt::tab_options( - heading.title.font.size = "1.5em", - heading.title.font.weight = "bolder", - column_labels.font.weight = "normal" - ) |> - gt::tab_style( - style = gt::cell_text(weight = "bold"), - locations = gt::cells_column_spanners() - ) |> - gt::tab_style( - style = gt::cell_text(align = "left"), - locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) - ) |> - gt::tab_style( - style = gt::cell_text(align = "right"), - locations = gt::cells_column_spanners(spanners = "outcome_span") - ) |> - gt::tab_footnote(footnote = pre_title_footnote_j, - locations = gt::cells_title()) |> - gt::tab_footnote(footnote = pre_units_footnote_j, - locations = gt::cells_column_spanners(spanners = "unit_span")) |> - gt::tab_footnote(footnote = pre_outcome_footnote_j, - locations = gt::cells_column_spanners(spanners = "outcome_span")) |> - gt::tab_footnote(footnote = pre_covariate_footnote_j, - locations = gt::cells_column_spanners(spanners = "covariates_span")) |> - gt::opt_css( - css = " - #preceptor_tbl_justice .gt_footnote { - max-width: 1px; - word-break: break-word; - } - " - ) |> - gt::as_raw_html(inline_css = FALSE) - -cat( - "```{=html}\n", - '
', - pre_gt_html_j, - "
\n", - "```\n", - sep = "" -) -``` - - -```{r} -#| echo: false -#| results: asis -data_tibble <- tibble::tribble( - ~`Respondent`, ~`Age`, ~`Height (cm)`, ~`Sex`, - "Emily Tran", "21", "162", "Female", - "James Okafor", "25", "178", "Male", - "Rachel Goldberg", "19", "159", "Female", - "...", "...", "...", "...", - "Tyler Brooks", "23", "170", "Male" -) - -data_title_footnote <- "National Health and Nutrition Examination Survey (NHANES), Centers for Disease Control and Prevention. Continuous survey data, 1999--present. We restrict to respondents aged 18 through 27 with non-missing height and sex --- about 1,300 rows. The original survey used a complex, stratified, multistage sampling design; we ignore the survey weights." -data_units_footnote <- "Each row is one NHANES participant examined between 1999 and the most recent release cycle. NHANES public files are anonymized, so the example rows use invented but plausible names." -data_outcome_footnote <- "Height in centimeters, measured by a trained examiner in NHANES's Mobile Examination Center using a fixed stadiometer. Each participant is measured once per survey visit with shoes removed." -data_covariate_footnote <- "Sex, recorded by the NHANES interviewer as 'Male' or 'Female'. NHANES does not separately record gender identity; the two-category coding dates from survey waves going back decades." - -data_gt_html <- gt::gt(data_tibble, id = "data_tbl") |> - gt::tab_header(title = "Data: NHANES Young Adults, Ages 18–27") |> - gt::tab_spanner(label = "Unit", id = "unit_span", - columns = c(`Respondent`, `Age`)) |> - gt::tab_spanner(label = "Outcome", id = "outcome_span", - columns = c(`Height (cm)`)) |> - gt::tab_spanner(label = "Covariate", id = "covariates_span", - columns = c(`Sex`)) |> - gt::cols_align(align = "left" , columns = c(`Respondent`, `Sex`)) |> - gt::cols_align(align = "right", columns = c(`Age`, `Height (cm)`)) |> - gt::fmt_markdown(columns = gt::everything()) |> - gt::tab_options( - heading.title.font.size = "1.5em", - heading.title.font.weight = "bolder", - column_labels.font.weight = "normal" - ) |> - gt::tab_style( - style = gt::cell_text(weight = "bold"), - locations = gt::cells_column_spanners() - ) |> - gt::tab_style( - style = gt::cell_text(align = "left"), - locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) - ) |> - gt::tab_style( - style = gt::cell_text(align = "right"), - locations = gt::cells_column_spanners(spanners = "outcome_span") - ) |> - gt::tab_footnote(footnote = data_title_footnote, - locations = gt::cells_title()) |> - gt::tab_footnote(footnote = data_units_footnote, - locations = gt::cells_column_spanners(spanners = "unit_span")) |> - gt::tab_footnote(footnote = data_outcome_footnote, - locations = gt::cells_column_spanners(spanners = "outcome_span")) |> - gt::tab_footnote(footnote = data_covariate_footnote, - locations = gt::cells_column_spanners(spanners = "covariates_span")) |> - gt::opt_css( - css = " - #data_tbl .gt_footnote { - max-width: 1px; - word-break: break-word; - } - " - ) |> - gt::as_raw_html(inline_css = FALSE) - -cat( - "```{=html}\n", - '
', - data_gt_html, - "
\n", - "```\n", - sep = "" -) -``` - -### Exercise 1 - -In your own words, describe the components of Justice when working on a data science problem. - -```{r justice-1} -question_text(NULL, - message = "Justice reviews the Population Table and selects the formula for the data generating mechanism.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Justice is about concerns that you (or your critics) might have, reasons why the model you create might not work as well as you hope. The exercises that follow walk through several such concerns by name. Justice also picks the formula for the data generating mechanism --- different outcome variables call for different formulas, a piece that gets more technical names in later tutorials. - -### Exercise 2 - -In your own words, define "validity" as we use the term. - -```{r justice-2} -question_text(NULL, - message = "Validity is the consistency, or lack thereof, in the columns of the data set and the corresponding columns in the Preceptor Table.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Validity is always about the columns in the Preceptor Table and the data. Just because columns from these two different tables have the same name does not mean that they are the same thing. - -### Exercise 3 - -Provide one reason why the assumption of validity might not hold for the outcome variable `height` or for one of the covariates. Use the words "column" or "columns" in your answer. - -```{r justice-3} -question_text(NULL, - message = "The `height` column in the data is recorded by NHANES's trained examiners using a fixed stadiometer, with shoes removed. The `height` column in the Preceptor Table is recorded at USMC enlistment with whatever instrument the recruiting office has on hand --- possibly with shoes on, possibly under different instructions. Even though both columns are labeled 'Height (cm),' the two measurement procedures are not the same, so the columns may not represent the same underlying quantity.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -In order to consider the Preceptor Table and the data to be drawn from the same population, the columns from one must have a *valid correspondence* with the columns in the other. Validity, if true (or at least reasonable), allows us to construct the Population Table, which is the first step in Justice. - -### Exercise 4 - -In your own words, define a Population Table. - -```{r justice-4} -question_text(NULL, - message = "The Population Table includes a row for each unit/time combination in the underlying population from which both the Preceptor Table and the data are drawn.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -The Population Table is almost always much bigger than the combination of the Preceptor Table and the data, because if we can really assume that both are part of the same population, then that population must cover a broad universe of time and units. - -### Exercise 5 - -Specify the unit/time combinations which define each row in this Population Table. - -```{r justice-5} -question_text(NULL, - message = "Each row is one young adult (aged 18--27) in a given year. Data rows are NHANES participants examined between 1999 and the most recent release. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Other rows cover all the young adults alive during the overlapping span who appear in neither table.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -```{r} -#| echo: false -#| results: asis -pop_tibble <- tibble::tribble( - ~Source , ~`Young Adult` , ~`Year`, ~`Height (cm)`, ~`Sex` , - "..." , "..." , "..." , "..." , "..." , - "Data" , "Emily Tran" , "2002" , "162" , "Female" , - "Data" , "James Okafor" , "2005" , "178" , "Male" , - "Data" , "..." , "..." , "..." , "..." , - "Data" , "Tyler Brooks" , "2024" , "170" , "Male" , - "..." , "..." , "..." , "..." , "..." , - "Preceptor", "Maya Castillo" , "2026" , "150" , "Female" , - "Preceptor", "Tyler Brooks" , "2026" , "170" , "Male" , - "Preceptor", "..." , "..." , "..." , "..." , - "Preceptor", "Grace Sullivan" , "2026" , "160" , "Female" , - "..." , "..." , "..." , "..." , "..." -) - -pop_title_footnote <- "This table combines NHANES adult-height data from 1999 onward with the Preceptor Table's 2026 USMC recruits, drawn from the same broader population of US young adults aged 18--27." -pop_units_footnote <- "Each row is one young adult at one point in time. Data rows are NHANES survey participants between the ages of 18 and 27. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Tyler Brooks appears in both blocks --- surveyed by NHANES in 2024 at age 23, enlisting in 2026 at 25. Same person, different times, different rows." -pop_outcome_footnote <- "Height in centimeters. NHANES measures height directly in the Mobile Examination Center; USMC records height at enlistment. The two measurement procedures are close enough that we treat them as the same column, subject to the validity concern raised in Exercise 3." -pop_covariate_footnote <- "Sex, recorded as Male or Female in both data sources." - -pop_gt_html <- gt::gt(pop_tibble, id = "population_tbl") |> - gt::tab_header(title = "Population Table") |> - gt::tab_spanner(label = "Unit/Time", id = "unit_span", - columns = c(`Young Adult`, `Year`)) |> - gt::tab_spanner(label = "Outcome" , id = "outcome_span", - columns = c(`Height (cm)`)) |> - gt::tab_spanner(label = "Covariate", id = "covariates_span", - columns = c(`Sex`)) |> - gt::cols_align(align = "left" , columns = c(`Source`, `Young Adult`, `Sex`)) |> - gt::cols_align(align = "right", columns = c(`Year`, `Height (cm)`)) |> - gt::fmt_markdown(columns = gt::everything()) |> - gt::tab_options( - heading.title.font.size = "1.5em", - heading.title.font.weight = "bolder", - column_labels.font.weight = "normal" - ) |> - gt::tab_style( - style = gt::cell_text(weight = "bold"), - locations = gt::cells_column_spanners() - ) |> - gt::tab_style( - style = gt::cell_text(align = "left"), - locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) - ) |> - gt::tab_style( - style = gt::cell_text(align = "right"), - locations = gt::cells_column_spanners(spanners = "outcome_span") - ) |> - gt::tab_footnote(footnote = pop_title_footnote, - locations = gt::cells_title()) |> - gt::tab_footnote(footnote = pop_units_footnote, - locations = gt::cells_column_spanners(spanners = "unit_span")) |> - gt::tab_footnote(footnote = pop_outcome_footnote, - locations = gt::cells_column_spanners(spanners = "outcome_span")) |> - gt::tab_footnote(footnote = pop_covariate_footnote, - locations = gt::cells_column_spanners(spanners = "covariates_span")) |> - gt::opt_css( - css = " - #population_tbl .gt_footnote { - max-width: 1px; - word-break: break-word; - } - " - ) |> - gt::as_raw_html(inline_css = FALSE) - -cat( - "```{=html}\n", - '
', - pop_gt_html, - "
\n", - "```\n", - sep = "" -) -``` - -### Exercise 6 - -In your own words, define the assumption of "stability" when employed in the context of data science. - -```{r justice-6} -question_text(NULL, - message = "Stability means that the relationship between the columns in the Population Table is the same for three categories of rows: the data, the Preceptor Table, and the larger population from which both are drawn.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -The most natural way for stability to fail is *time*. Our data was gathered in earlier NHANES cycles, but our recruits will enlist next year. Stability holds only if the world that produced the data behaves the same as the world the recruits enter. That's not guaranteed here: NHANES spans 25 years, and compared to next year's recruit cohort, that's a long enough window for the relationship between sex and adult height to drift. -We will sharpen this idea --- the *parameters* of the model, not the levels of any single variable --- in the next exercise. - -### Exercise 7 - -Provide one reason why the assumption of stability might not be true in this case. - -```{r justice-7} -question_text(NULL, - message = "The *intercept* of the height-on-sex regression --- the expected height of a female young adult --- has likely drifted upward over the 25-year span of NHANES, while the *slope* (the male-female height gap, the `sexMale` coefficient) has stayed roughly constant. By 2026 the underlying $\\beta_0$ for our recruits is probably a centimeter or two higher than the average across NHANES's reporting years. That is a *parameter* shift; if it is large enough, stability is violated. The fact that average heights in the population have changed is not by itself the issue --- a distribution change in `height` is something a stable model can absorb. What stability requires is that $\\beta_0$ and $\\beta_1$ stay put.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -A common confusion is to point at any change between two eras and call it a stability violation. It isn't. The distribution of `height` in the population may have shifted upward; the spread of heights may have narrowed; the mix of ethnicities in the recruiting pool may differ from NHANES participants. None of those, on their own, is a stability violation. What stability requires is that the *parameters* of our model --- $\beta_0$ and $\beta_1$ --- be the same in the data, in the Preceptor Table, and in the broader population. Distribution shifts are everywhere; parameter shifts are what hurt us. - -### Exercise 8 - -We use our data to make inferences about the overall population. We use information about the population to make inferences about the Preceptor Table: Data → Population → Preceptor Table. In your own words, define the assumption of "representativeness." - -```{r justice-8} -question_text(NULL, - message = "Representativeness, or the lack thereof, concerns two relationships among the rows in the Population Table. The first is between the data and the other rows. The second is between the other rows and the Preceptor Table.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Ideally, we would like both the Preceptor Table and our data to be random samples from the population. Sadly, this is almost never the case. - -### Exercise 9 - -We do not use the data directly to estimate missing values in the Preceptor Table. Instead, we use the data to learn about the overall population. Provide one reason, involving the relationship between the data and the population, why the assumption of representativeness might not be true in this case. - -```{r justice-9} -question_text(NULL, - message = "NHANES is a voluntary survey. People who agree to participate (and who show up for the Mobile Examination Center visit) are not a random sample of US young adults --- they skew toward people with time, transportation, and some trust in government research. On top of that, the `recruits` tibble we are fitting on is itself a non-random subsample of NHANES: 40 male rows and 10 female rows, deliberately unbalanced. Both layers --- the survey participation and our 40/10 choice --- mean the data we model is not a random sample of US young adults.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -When representativeness is violated, our parameter estimates *might* be biased --- not "will be," because by pure luck a non-representative sample can still produce the right number, but we have no principled reason to expect that to happen. There is a *second* representativeness relationship --- between the population and the Preceptor Table --- that later tutorials will exercise; for now, the data ↔ population link is what we focus on. - -### Exercise 10 - -A statistical model consists of two parts: the probability family and the link function. The probability family is the probability distribution that generates the randomness in our data. Our outcome variable is `height` --- continuous. What probability family should we use? - -```{r justice-10} -question_text(NULL, - message = "Normal: $Y \\sim N(\\mu, \\sigma^2)$. The Normal distribution is the standard choice for continuous outcomes.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 4) -``` - -### - -While you're here, also add `library(tidymodels)` to the setup chunk in `analysis.qmd` so the rest of Courage can call `linear_reg()` and `fit()`. - -### Exercise 11 - -The link function is the mathematical formula that links our data to the unknown parameters. What is the link function for a continuous outcome? - -```{r justice-11} -question_text(NULL, - message = "Linear: $\\mu = \\beta_0 + \\beta_1 X_1 + \\ldots + \\beta_n X_n$. Our model will be a linear regression.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 4) -``` - -### - -While you're here, also add `library(broom)` to the setup chunk in `analysis.qmd` so we can call `tidy()` on the fitted model in Courage. - -### Exercise 12 - -Write one sentence highlighting a potential weakness in your model. Derive it from possible problems with the assumptions above. We will add this to our summary paragraph. So far our version of the summary paragraph looks like this: - -> People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year. - -Your version will be somewhat different. - -```{r justice-12} -question_text(NULL, - message = "One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Add a weakness sentence to the summary paragraph in your QMD. You can modify your paragraph, but don't copy/paste our answer exactly. Render the QMD, then commit/push. - - -## Courage -### - -Justice gives us the Population Table and the abstract data generating mechanism. - -```{r} -#| echo: false -#| results: asis -pop_tibble_c <- tibble::tribble( - ~Source , ~`Young Adult` , ~`Year`, ~`Height (cm)`, ~`Sex` , - "..." , "..." , "..." , "..." , "..." , - "Data" , "Emily Tran" , "2002" , "162" , "Female" , - "Data" , "James Okafor" , "2005" , "178" , "Male" , - "Data" , "..." , "..." , "..." , "..." , - "Data" , "Tyler Brooks" , "2024" , "170" , "Male" , - "..." , "..." , "..." , "..." , "..." , - "Preceptor", "Maya Castillo" , "2026" , "150" , "Female" , - "Preceptor", "Tyler Brooks" , "2026" , "170" , "Male" , - "Preceptor", "..." , "..." , "..." , "..." , - "Preceptor", "Grace Sullivan" , "2026" , "160" , "Female" , - "..." , "..." , "..." , "..." , "..." -) - -pop_title_footnote_c <- "This table combines NHANES adult-height data from 1999 onward with the Preceptor Table's 2026 USMC recruits, drawn from the same broader population of US young adults aged 18--27." -pop_units_footnote_c <- "Each row is one young adult at one point in time. Data rows are NHANES survey participants between the ages of 18 and 27. Preceptor rows are the ~5,000 USMC recruits expected to enlist in the coming year. Tyler Brooks appears in both blocks." -pop_outcome_footnote_c <- "Height in centimeters. NHANES measures height directly in the Mobile Examination Center; USMC records height at enlistment." -pop_covariate_footnote_c <- "Sex, recorded as Male or Female in both data sources." - -pop_gt_html_c <- gt::gt(pop_tibble_c, id = "population_tbl_courage") |> - gt::tab_header(title = "Population Table") |> - gt::tab_spanner(label = "Unit/Time", id = "unit_span", - columns = c(`Young Adult`, `Year`)) |> - gt::tab_spanner(label = "Outcome" , id = "outcome_span", - columns = c(`Height (cm)`)) |> - gt::tab_spanner(label = "Covariate", id = "covariates_span", - columns = c(`Sex`)) |> - gt::cols_align(align = "left" , columns = c(`Source`, `Young Adult`, `Sex`)) |> - gt::cols_align(align = "right", columns = c(`Year`, `Height (cm)`)) |> - gt::fmt_markdown(columns = gt::everything()) |> - gt::tab_options( - heading.title.font.size = "1.5em", - heading.title.font.weight = "bolder", - column_labels.font.weight = "normal" - ) |> - gt::tab_style( - style = gt::cell_text(weight = "bold"), - locations = gt::cells_column_spanners() - ) |> - gt::tab_style( - style = gt::cell_text(align = "left"), - locations = gt::cells_column_spanners(spanners = c("unit_span", "covariates_span")) - ) |> - gt::tab_style( - style = gt::cell_text(align = "right"), - locations = gt::cells_column_spanners(spanners = "outcome_span") - ) |> - gt::tab_footnote(footnote = pop_title_footnote_c, - locations = gt::cells_title()) |> - gt::tab_footnote(footnote = pop_units_footnote_c, - locations = gt::cells_column_spanners(spanners = "unit_span")) |> - gt::tab_footnote(footnote = pop_outcome_footnote_c, - locations = gt::cells_column_spanners(spanners = "outcome_span")) |> - gt::tab_footnote(footnote = pop_covariate_footnote_c, - locations = gt::cells_column_spanners(spanners = "covariates_span")) |> - gt::opt_css( - css = " - #population_tbl_courage .gt_footnote { - max-width: 1px; - word-break: break-word; - } - " - ) |> - gt::as_raw_html(inline_css = FALSE) - -cat( - "```{=html}\n", - '
', - pop_gt_html_c, - "
\n", - "```\n", - sep = "" -) -``` - -$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_n X_n + \epsilon$$ - -with $\epsilon \sim N(0, \sigma^2)$. - -### Exercise 1 - -In your own words, describe the components of the virtue of Courage for analyzing data. - -```{r courage-1} -question_text(NULL, - message = "Courage creates the data generating mechanism.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -The abstract data generating mechanism shown in the preamble uses generic variables --- $Y$, $X_1$, $X_2$, and so on --- because our purpose is to describe the general mathematical structure. At the end of Justice we have decided the functional family (linear with normal errors), but not the number or identity of covariates we will ultimately include. Courage fills in the missing pieces. - -### Exercise 2 - -Because our outcome variable is continuous, add a code chunk to `analysis.qmd` that begins building the model with `linear_reg(engine = "lm")`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r courage-2-test} -#| echo: true -linear_reg(engine = "lm") -``` - -### - -**[tidymodels](https://www.tidymodels.org/)** is the most popular framework for estimating models among R users. *[Tidy Modeling with R](https://www.tmwr.org/)* by Max Kuhn and Julia Silge is a great introduction. - -### Exercise 3 - -Extend the chunk in `analysis.qmd` to `linear_reg(engine = "lm") |> fit(height ~ sex, data = recruits)`. This is our first candidate model: outcome on the binary `sex`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-3} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r courage-3-test} -#| echo: true -linear_reg(engine = "lm") |> - fit(height ~ sex, data = recruits) -``` - -### - -A two-level character variable like `sex` is turned into a 0/1 "dummy" variable, renamed `sexMale` (the level for which the dummy is 1). Math doesn't speak text; dummies do the translation. - -### Exercise 4 - -Here is the parameter table for this candidate model: - -```{r} -#| echo: false -linear_reg() |> - set_engine("lm") |> - fit(height ~ sex, data = recruits) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 1))) -``` - -Write a sentence interpreting the `sexMale` estimate of 8.9. - -```{r courage-4} -question_text(NULL, - message = "When we compare two groups of recruits differing only in sex, the male group has an expected height about 8.9 centimeters taller than the female group. The intercept of 165.9 is the expected height of a female recruit; adding `sexMale` (= 1 for male) gives the expected male height of about 174.8 cm.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -Notice the language: *"when we compare two groups,"* not *"raising sex by one unit increases height."* For a binary covariate the comparison framing is the only honest reading. The same comparison framing applies to continuous covariates too --- which we can see by trying one. The 95% confidence interval for `sexMale` is [4.3, 13.4] centimeters; zero is well outside that interval, so the model says the *expected* height of male recruits is taller than the expected height of female recruits. That does not mean every male is taller than every female --- the two distributions overlap, and there are plenty of female recruits taller than some male recruits. - -### Exercise 5 - -Change the formula to `fit(height ~ age, data = recruits)`. This refits with `age`, a continuous covariate, instead of `sex`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r courage-5-test} -#| echo: true -linear_reg(engine = "lm") |> - fit(height ~ age, data = recruits) -``` - -### - -Continuous covariates produce a single slope coefficient. The reading is per-unit comparison: how much do two groups, differing by one unit of the covariate, differ in expected outcome? - -### Exercise 6 - -Here is the parameter table for the second candidate model: - -```{r} -#| echo: false -linear_reg() |> - set_engine("lm") |> - fit(height ~ age, data = recruits) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 1))) -``` - -Write a sentence interpreting the slope on `age`. - -```{r courage-6} -question_text(NULL, - message = "The slope on `age` is essentially zero (0.0 cm per year, with a 95% confidence interval that crosses zero in both directions). When we compare two groups of recruits differing by one year of age (within our 18--27 window), their expected heights are indistinguishable. Adult height has stopped changing by the late teens; this fit confirms that.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -The wide confidence interval that includes zero is exactly what *"this covariate doesn't help"* looks like in a parameter table. We tried `age`, and it brought no signal beyond what `sex` already gives us. Our final model is therefore `height ~ sex` --- the candidate from Exercise 4 --- which is the model the rest of the tutorial uses. - -### Exercise 7 - -We've settled on the final fit. In `analysis.qmd`, add a new code chunk that assigns the result of `linear_reg() |> set_engine("lm") |> fit(height ~ sex, data = recruits)` to `fit_recruits`, then on the next line prints `fit_recruits`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-7} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r courage-7-test} -#| echo: true -fit_recruits <- linear_reg() |> - set_engine("lm") |> - fit(height ~ sex, data = recruits) -fit_recruits -``` - -### - -In data science, **we deal with words, math, and code, but the most important of these is code.** We created the mathematical structure of the model and then wrote a model formula in order to estimate the unknown parameters. - -### Exercise 8 - -Re-fitting the model on every render is wasteful --- the fit's the same each time. `fit_recruits` earns its own cached create chunk. The next three exercises set that up: clean the chunk, turn caching on, add the cache directory to `.gitignore`. - -First, clean up. In your chunk from Exercise 7, delete the `fit_recruits` line at the bottom, leaving just the `fit_recruits <- linear_reg() ...` assignment. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-8} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -
fit_recruits <- linear_reg() |>
-  set_engine("lm") |>
-  fit(height ~ sex, data = recruits)
-
- -### - -A chunk ready to cache should contain only the artifact it builds --- nothing else. - -### Exercise 9 - -Add `#| cache: true` as the first line inside the chunk. Render the QMD. Rendering creates an `analysis_cache/` directory next to your `analysis.qmd`, which holds the saved fit so Quarto can reload it on subsequent renders rather than re-fitting every time. To confirm the directory is there, run `ls` in the **bash Terminal**. CP/CR. - -```{r courage-9} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -
$ ls
-analysis.qmd  analysis_cache  analysis_files  analysis.html
-
- -### - -For a 50-row linear regression the savings are small; for the larger model fits in later tutorials, caching can be the difference between a one-second and a one-minute render. - -### Exercise 10 - -The Source Control change count just jumped, because `analysis_cache/` and everything inside it now look like untracked work to git. Cached files should never go to GitHub --- they're machine-specific and regenerated on every render. - -Add `analysis_cache` to your `.gitignore` on its own line. In the R Terminal, run `show_file(".gitignore")`. CP/CR. - -```{r courage-10} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -
analysis_files
-analysis_cache
-
- -### - -The Source Control change count should drop back to whatever it was before the render. Cached objects live next to the QMD, not in version control. - -With the parameters estimated, we can write the fitted model as a concrete formula: - -$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ - -with residuals drawn from $N(0, 41.0)$ --- a residual standard deviation of about 6.4 cm. - -Three differences from the abstract form shown in the preamble. First, the parameters are replaced by their best estimates (rounded to one decimal place, fine enough that the estimates and their confidence-interval bounds remain visibly distinct). Second, the error term is gone --- this version of the DGM answers an expected-value question, so the randomness sits to the side in the residuals line until a question needs it. Third, the left-hand side has a hat, because hats mark estimated values. - -**This is our data generating mechanism.** A data generating mechanism is just a formula, something which we can write down and implement with computer code. The randomness built into the DGM lives in the residuals line. - -### Exercise 11 - -Start a new working chunk below the cached `fit_recruits` chunk that runs `tidy(fit_recruits, conf.int = TRUE)`. This returns 95% intervals for all the parameters in our model. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r courage-11} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r courage-11-test} -#| echo: true -tidy(fit_recruits, conf.int = TRUE) -``` - -### - -`tidy()` is part of the **[broom](https://broom.tidymodels.org/)** package, used to summarize information from a wide variety of models. - -Here is a nicer-looking version of that parameter table, rendered with **[knitr](https://yihui.org/knitr/)**: - -```{r} -#| echo: false -library(knitr) - -tidy(fit_recruits, conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 1))) |> - kable( - caption = "Estimates and 95% Confidence Intervals for Parameters. Source: 50-row sample from NHANES via primer.tutorials::recruits" - ) -``` - -Rounding to one decimal place keeps the displayed values consistent with the fitted DGM formula above and is fine enough that the estimate and the confidence-interval bounds are visibly distinguishable in every row. A blanket "three significant figures" would risk collapsing values close to one another --- like an intercept of 165.9 and its upper bound 170.0 --- into ambiguous displays. You need rounding finer than the standard error to keep adjacent values apart. At the very least, a table like this should include a title and a caption with the data source. - -### - -The more variables we add, the more difficult it is to interpret the meaning of any particular coefficient. But interpretation also becomes less important. We don't really care about coefficients. We care about using our model to estimate quantities of interest. - -### Exercise 12 - -Add a sentence to your project summary explaining the structure of the model. Something like: "We model XX [concept of outcome, not variable name], [values of XX], as a [linear/logistic/multinomial/ordinal] function of XX [and maybe other covariates]." - -Recall the beginning of our version of the summary: - -> People vary in height, and that variation is strongly patterned by sex. Using data from the National Health and Nutrition Examination Survey (NHANES) conducted by the Centers for Disease Control and Prevention, we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. - -```{r courage-12} -question_text(NULL, - message = "We model height as a normally distributed variable which is a linear function of sex.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Read our answer. Do not copy/paste exactly. Add your two sentences to the summary paragraph. Render the QMD. - - - -## Temperance -### - -Courage provides the data generating mechanism. We can express the DGM in four ways: - -**In words:** We model height as a normally distributed variable which is a linear function of sex. - -**In R code:** - -```r -fit_recruits <- linear_reg() |> - set_engine("lm") |> - fit(height ~ sex, data = recruits) -``` - -**In a parameter table:** - -```{r} -#| echo: false -fit_recruits |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 1))) -``` - -**As a mathematical formula:** - -$$\widehat{\text{height}} = 165.9 + 8.9 \cdot \text{sexMale}$$ - -with residuals drawn from $N(0, 41.0)$ --- a residual standard deviation of about 6.4 cm. The hatted equation is the expected-value mode of the DGM --- plug in a sex, get an expected height --- while the residuals line keeps the draw mode's randomness at hand for questions about individual recruits. - -### Exercise 1 - -In your own words, describe the use of Temperance in data science. - -```{r temperance-1} -question_text(NULL, - message = "Temperance interprets the data generating mechanism and then uses it to answer, with the help of graphics, the question(s) with which we began. Humility reminds us that this answer is always false.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Courage gave us the data generating mechanism, including the parameter values that Temperance now puts to work. The tool for that work is the **[marginaleffects](https://marginaleffects.com/)** package, with companion book *[Model to Meaning](https://marginaleffects.com/)* by Vincent Arel-Bundock. - -### Exercise 2 - -In the end, we don't really care about parameters, much less how to interpret them. Parameters are *imaginary*, like unicorns. We care about answers to our questions. *In the modern world, all parameters are nuisance parameters.* - -Add `library(marginaleffects)` to the setup chunk in `analysis.qmd` and render. CP/CR. - -```{r temperance-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - - -### - -

-#| message: false
-library(tidyverse)
-library(primer.tutorials)
-library(tidymodels)
-library(broom)
-library(marginaleffects)
-
- -### - -We should be modest in the claims we make. The posteriors we create are never the "truth." The assumptions we made to create the model are never perfect. Yet decisions made with flawed posteriors are almost always better than decisions made without them. - -### Exercise 3 - -What is the specific question we are trying to answer? - -```{r temperance-3} -question_text(NULL, - message = "What is the expected height of male and female USMC recruits next year?", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 3) -``` - -### - -Data science projects begin with a decision which we face. To make that decision wisely, we would like to have good estimates of many unknown numbers. Yet, in order to make progress, we need to drill down to one specific question. This leads to the creation of a data generating mechanism, which can then be used to answer lots of questions. - -### Exercise 4 - -In `analysis.qmd`, add a chunk that runs `predictions(fit_recruits)`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r temperance-4} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r temperance-4-test} -#| echo: true -predictions(fit_recruits) -``` - -### - -`predictions()` returns a data frame with one row for each observation in the data set used to fit the model. The 50 rows correspond to our 50-row `recruits` sample. Each row's `estimate` is that row's model-predicted height. Note already that the standard error column has two distinct values --- about 1.01 for male rows and about 2.02 for female rows. We will come back to this. - -### - -One DGM answers many questions. Ask it for an expected value and it is a deterministic formula: same sex in, same number out --- which is why these 50 rows show only two distinct estimates. Ask it what will happen to a single recruit and the randomness matters: individual heights have variance baked in, so the DGM must draw from its residual distribution, and two recruits of the same sex can end up with very different heights. The question, not habit, picks the mode. - -### Exercise 5 - -In `analysis.qmd`, add a chunk that runs `plot_predictions(fit_recruits, condition = "sex")`. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r temperance-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r temperance-5-test} -#| echo: true -plot_predictions(fit_recruits, condition = "sex") -``` - -### - -You should see two point estimates (one for Female, one for Male) with 95% confidence intervals. The gap between them is the model's estimate of the average height difference by sex. The Female interval is visibly wider than the Male interval --- the predicted height for a female recruit is less certain than the predicted height for a male recruit. That is the consequence of having drawn 10 female rows but 40 male rows: the standard error of a group mean shrinks roughly with the square root of the group's sample size, so quadrupling the sample size halves the standard error, and the female interval ends up about twice as wide as the male interval. - -### Exercise 6 - -Modify the previous chunk to add `draw = FALSE` to the `plot_predictions()` call. Render. In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r temperance-6} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r temperance-6-test} -#| echo: true -plot_predictions(fit_recruits, condition = "sex", draw = FALSE) -``` - -### - -Because `plot_predictions()` returns a ggplot object, you can continue with ggplot commands like `labs()`. But it can be useful to see the underlying values in the tibble and build your own plot directly. - -### Exercise 7 - -Create a beautiful plot starting from the output of `plot_predictions(fit_recruits, condition = "sex", draw = FALSE)`. Do this in your QMD (much easier than typing in the R Terminal directly). Title: the key variables. Subtitle: the key takeaway. Caption: data source. Axis labels: human-readable with units. - -Paste the plot code below. - -```{r temperance-7} -question_text(NULL, - message = " plot_predictions(fit_recruits, condition = \"sex\", draw = FALSE) |>\n ggplot(aes(x = sex, y = estimate)) +\n geom_pointrange(aes(ymin = conf.low, ymax = conf.high),\n color = \"grey20\", size = 1.1) +\n labs(\n title = \"Forecast Height for USMC Recruits\",\n subtitle = \"Male recruits are expected to be about 9 cm taller\",\n x = \"Sex\",\n y = \"Expected Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 12) -``` - -### - -Here is our version: - -```{r} -#| echo: false -#| message: false -#| warning: false -plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |> - ggplot(aes(x = sex, y = estimate)) + - geom_pointrange(aes(ymin = conf.low, ymax = conf.high), - color = "grey20", size = 1.1) + - labs( - title = "Forecast Height for USMC Recruits", - subtitle = "Male recruits are expected to be about 9 cm taller than female recruits", - x = "Sex", - y = "Expected Height (cm)", - caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" - ) + - theme_minimal() -``` - -Data science often involves this back-and-forth style of work. First, make a single chunk of code --- say, a new plot --- work well. This requires interactive work between the QMD and the R Terminal. Second, ensure that the entire QMD runs correctly on its own. - -### Exercise 8 - -Finalize the new graphics chunk in your QMD. Render to ensure it all works. In the R Terminal, run: - -``` -show_file("analysis.qmd", chunk = "Last") -``` - -CP/CR. - -```{r temperance-8} -question_text(NULL, - message = " > show_file(\"analysis.qmd\", chunk = \"Last\")\n ```{r}\n plot_predictions(fit_recruits, condition = \"sex\", draw = FALSE) |>\n ggplot(aes(x = sex, y = estimate)) +\n geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +\n labs(\n title = \"Forecast Height for USMC Recruits\",\n subtitle = \"Male recruits are expected to be about 9 cm taller\",\n x = \"Sex\", y = \"Expected Height (cm)\",\n caption = \"Source: 50-row sample from NHANES via primer.tutorials::recruits\"\n ) +\n theme_minimal()\n ```", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 8) -``` - -### - -Always remember: the map is not the territory. A beautiful graphic tells a story, but that story is always an imperfect representation of reality. Our models depend on assumptions that are never completely true. - -### Exercise 9 - -Write the last sentence of your summary paragraph. It describes at least one Quantity of Interest and a measure of uncertainty. It is OK if this QoI differs from the one you began with. - -```{r temperance-9} -question_text(NULL, - message = "We estimate that the average male USMC recruit will be about 175 cm tall and the average female recruit about 166 cm, with the female interval visibly wider than the male interval because we have far fewer female observations in our sample.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 5) -``` - -### - -Add a final sentence to your summary paragraph, but don't copy/paste our answer exactly. Render the QMD. - -### Exercise 10 - -Write a few sentences explaining why the estimates for the quantities of interest, and the uncertainty, might be wrong. Suggest alternative estimates and a wider confidence interval if warranted. - -```{r temperance-10} -question_text(NULL, - message = "USMC recruits are more physically fit and probably taller than the general young-adult population measured by NHANES. A better point estimate would nudge male recruit heights upward by a centimeter or two and female recruit heights similarly. The confidence intervals should be wider than what the model reports --- our 50-row sample is small (only 10 of those rows are female), and the bigger sources of uncertainty (the NHANES-to-USMC validity gap and the voluntary-participation representativeness gap) aren't captured by the model's standard error.", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -Always go back to your Preceptor Table --- the information which, if you had it, would make answering your question easy. In almost all real-world cases, the Preceptor Table and the data are fairly different. So, even a perfectly estimated statistical model is rarely as useful as we might like. - -### Exercise 11 - -Rearrange the material in your QMD so the order is graphic, then summary paragraph. The chunk that creates `fit_recruits` must occur before the chunk that creates the graphic. You can keep or discard the math at your discretion. Render the QMD. In the R Terminal, run: - -``` -show_file("analysis.qmd") -``` - -CP/CR. - -```{r temperance-11} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 10) -``` - -### - -

----
-title: "Recruits"
-author: "<your name>"
-execute:
-  echo: false
----
-
-#| message: false
-library(tidyverse)
-library(tidymodels)
-library(primer.tutorials)
-library(broom)
-library(marginaleffects)
-
-#| cache: true
-fit_recruits <- linear_reg() |>
-  set_engine("lm") |>
-  fit(height ~ sex, data = recruits)
-
-plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |>
-  ggplot(aes(x = sex, y = estimate)) +
-  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
-  labs(
-    title = "Forecast Height for USMC Recruits",
-    subtitle = "Male recruits are expected to be about 9 cm taller",
-    x = "Sex", y = "Expected Height (cm)",
-    caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits"
-  ) +
-  theme_minimal()
-
-People vary in height, and that variation is strongly patterned by sex. Using a 50-row teaching sample drawn from the National Health and Nutrition Examination Survey (NHANES, conducted by the Centers for Disease Control and Prevention), we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. We model height as a normally distributed variable which is a linear function of sex. We estimate that the average male recruit will be about 175 cm tall and the average female recruit about 166 cm, with confidence intervals that are visibly wider for the female estimate because there are only ten female rows in our sample.
-
- -### - -This is the version of your QMD file at which your teacher is most likely to look closely. - -### Exercise 12 - -Publish your rendered QMD to GitHub Pages. In the bash Terminal (not the R Terminal!), run: - -``` -quarto publish gh-pages analysis.qmd -``` - -Copy/paste the resulting URL below. - -```{r temperance-12} -question_text(NULL, - message = "https://.github.io/recruits/analysis.html", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 2) -``` - -### - -Commit/push everything. - -### Exercise 13 - -Copy/paste the URL to your GitHub repo. - -```{r temperance-13} -question_text(NULL, - message = "https://github.com//recruits", - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 2) -``` - -### - -We can never know all the entries in the Preceptor Table. That knowledge is reserved for God. If all our assumptions are correct, then our DGM is true --- it accurately describes the way in which the world works. There is no better way to predict the future, or to model the past, than to use it. Sadly, this will only be the case with toy examples involving things like coins and dice. - -The world confronts us. Make decisions we must. - - -## Summary -### - -```{r} -#| echo: false -#| message: false -#| warning: false -plot_predictions(fit_recruits, condition = "sex", draw = FALSE) |> - ggplot(aes(x = sex, y = estimate)) + - geom_pointrange(aes(ymin = conf.low, ymax = conf.high), - color = "grey20", size = 1.1) + - labs( - title = "Forecast Height for USMC Recruits", - subtitle = "Male recruits are expected to be about 9 cm taller than female recruits", - x = "Sex", - y = "Expected Height (cm)", - caption = "Source: 50-row sample from NHANES via primer.tutorials::recruits" - ) + - theme_minimal() -``` - -People vary in height, and that variation is strongly patterned by sex. Using a 50-row teaching sample drawn from the National Health and Nutrition Examination Survey (NHANES, conducted by the Centers for Disease Control and Prevention), we seek to estimate the average height of male and female USMC recruits for the coming year. One concern is that USMC recruits are more physically fit than the general young-adult population and so may be systematically taller, meaning that the NHANES-based estimates probably understate recruit heights. We model height as a normally distributed variable which is a linear function of sex. We estimate that the average male recruit will be about 175 cm tall and the average female recruit about 166 cm, with the female interval visibly wider than the male interval because we have only ten female rows in our sample. - -This forecast would help you estimate how many small-versus-large uniforms to order. It does not address every question the logistics office will face --- it says nothing about weight, shoulder width, or inseam. But it gives a defensible starting point. - -The world is always more uncertain than our models would have us believe. - -```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} -``` diff --git a/inst/tutorials/dgm-linear/tutorial.Rmd b/inst/tutorials/dgm-linear/tutorial.Rmd deleted file mode 100644 index 621f911..0000000 --- a/inst/tutorials/dgm-linear/tutorial.Rmd +++ /dev/null @@ -1,1257 +0,0 @@ ---- -title: DGM Linear -author: Var Kurapati -tutorial: - id: dgm-linear - version: 1.2.0 -output: - learnr::tutorial: - df_print: default - progressive: true - allow_skip: true -runtime: shiny_prerendered -description: 'Learn to fit and interpret linear regression models using tidymodels and marginaleffects, applied to housing prices in Ames, Iowa.' ---- - -```{r setup, include=FALSE} -library(learnr) -library(tutorial.helpers) -library(knitr) -library(tidyverse) -library(tidymodels) -library(marginaleffects) - -knitr::opts_chunk$set(echo = FALSE) -knitr::opts_chunk$set(out.width = '90%') -options(tutorial.exercise.timelimit = 600, - tutorial.storage = "local") - -mod_linear <- linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames) -``` - -```{r info-section, child = system.file("child_documents/info_section.Rmd", package = "tutorial.helpers")} -``` - -## Introduction -### - -The `ames` dataset records 2,930 houses sold in Ames, Iowa between 2006 and 2010 — the size, age, and condition of each house, along with the price it actually sold for. - -This tutorial asks a simple question: **can we predict a house's sale price from its size and age?** - -`Sale_Price` is a number — not a category — and that is information a model can use. We use two packages to build and interpret such a model. **[tidymodels](https://www.tidymodels.org/)** fits it. **[marginaleffects](https://marginaleffects.com/)** interprets it — turning raw coefficients into predicted dollar amounts that are immediately understandable. - -Interpreting a model with **marginaleffects** almost always comes down to one distinction: **marginal** vs. **conditional**. A marginal prediction runs the model on every *real* house in the dataset and averages the results (the `by` argument). A conditional prediction builds one or more *hypothetical* houses, holding some variables fixed and varying others (the `condition` argument). So every time you see `by`, you are averaging over real houses, and every time you see `condition`, you are looking at a hypothetical. - -Every dataset is a shadow cast by some real process — a **data generating mechanism**, or DGM — never the thing itself. Two moves lie ahead: *Courage*, proposing a structure and fitting it to real data (Tidymodels), and *Temperance*, using that fitted machine — with humility — to answer the question we started with (Marginal Effects). - -### Exercise 1 - -You should be connected to a repo named `dgm-linear`. If you are not, create one and connect to it. - -Create a new file, `analysis.qmd`, with the title `"Ames Housing"` and your name as the author. Add `execute: echo: false` to the YAML header. In a bash Terminal, render it: - -``` -quarto render analysis.qmd -``` - -Open `analysis.html` with Live Server. Create a `.gitignore` with `analysis_files` on the first line. Commit and push. - -In the R Terminal, run: - -``` -show_file(".gitignore") -``` - -CP/CR. - -```{r introduction-1} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 3) -``` - -### - -``` -analysis_files -``` - -### - -Every dataset comes from somewhere. `ames` is a record of houses that were actually sold in Ames, Iowa between 2006 and 2010, compiled by the local assessor's office — not a random sample of all houses everywhere. - -### Exercise 2 - -In a new code chunk in `analysis.qmd`, load the libraries you'll need. Add `#| message: false`. Render. - -``` -#| message: false -library(tidyverse) -library(tidymodels) -library(marginaleffects) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r introduction-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 10) -``` - -### - -
#| message: false
-library(tidyverse)
-library(tidymodels)
-library(marginaleffects)
-
- -### - -`ames` needs no loading of its own — it comes with **modeldata**, which `library(tidymodels)` attaches. The data is a shadow of the Ames housing market; the DGM we build will be our best account of the process that cast it. - -### Exercise 3 - -In a new code chunk, print the first 10 rows of `ames`, narrowed to the three columns this tutorial uses. Render. - -``` -ames |> - select(Sale_Price, Gr_Liv_Area, Year_Built) |> - head(10) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r introduction-3} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 4) -``` - -### - -```{r introduction-3-test} -#| echo: true -ames |> - select(Sale_Price, Gr_Liv_Area, Year_Built) |> - head(10) -``` - -### - -`Sale_Price` is the outcome we want to predict — the actual price each house sold for in dollars. `Gr_Liv_Area` is above-ground living area in square feet. `Year_Built` is the year the house was constructed. - -### Exercise 4 - -Ask AI to create a scatterplot of `Sale_Price` vs. `Gr_Liv_Area`. Paste its code into a new chunk. Render. - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r introduction-4} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 10) -``` - -### - -Here's our version, for comparison — yours doesn't need to match exactly, just the same idea: - -```{r introduction-4-test} -#| echo: true -ames |> - ggplot(aes(x = Gr_Liv_Area, y = Sale_Price)) + - geom_point(alpha = 0.3) + - labs( - title = "House sale price vs above-ground living area", - x = "Above-ground living area (sq ft)", - y = "Sale price ($)" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -### - -There is a clear positive relationship — larger houses sell for more. The relationship is roughly linear, which means a linear regression model is a reasonable choice for this data. - -### Exercise 5 - -Commit and push, with a message like "Add Ames housing data and initial exploration." In the bash Terminal, run: - -``` -git log -1 -``` - -CP/CR. - -```{r introduction-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -``` -commit a1b2c3d4e5f6... -Author: Var Kurapati -Date: Mon Jul 27 2026 - Add Ames housing data and initial exploration -``` - -### - -Before we fit anything, it's worth asking where this data came from. Whatever `mod_linear` learns, it learns about *this* time and place — Ames, Iowa, 2006–2010 — not housing markets in general. Keep asking, as we build the model ahead: *How did this data come to be?* - -## Tidymodels -### - -Every dataset is produced by some real-world process — a **data generating mechanism**, or DGM. Building one means committing to a structure, choosing a method for estimating it, and fitting it to data: **Structure → Method → Fitting**. `Sale_Price` is a number, so the structure here is `linear_reg()`. - -There is no single "correct" DGM for a given question. Which right-hand-side variables to include is a judgment call — one we will make, and remake, several times in this section. - -### Exercise 1 - -In a new code chunk, count how many houses fall into each neighborhood, sorted from most to least common. Render. - -``` -ames |> - count(Neighborhood, sort = TRUE) |> - head(10) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-1} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 6) -``` - -### - -```{r tidymodels-1-test} -#| echo: true -ames |> - count(Neighborhood, sort = TRUE) |> - head(10) -``` - -### - -North Ames is the most common neighborhood with 443 houses. The uneven distribution across neighborhoods is worth remembering — any model we fit will be better calibrated for neighborhoods with many houses than for those with few. - -### Exercise 2 - -In a new code chunk, specify a linear regression with `linear_reg()` and pipe it into `set_engine("lm")`. Print the result. Render. - -``` -linear_reg() |> - set_engine("lm") -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 4) -``` - -### - -```{r tidymodels-2-test} -#| echo: true -linear_reg() |> - set_engine("lm") -``` - -### - -`linear_reg()` proposes the **structure**: a numeric outcome predicted by a linear combination of predictors. `set_engine("lm")` picks the **method**: R's classical least-squares fitting. These are separate decisions — the DGM and our procedure for estimating it are separate things. - -### Exercise 3 - -Update your chunk to pipe into `fit(Sale_Price ~ Gr_Liv_Area, data = ames)` — just one predictor for now. Then pipe the result into `tidy(conf.int = TRUE)` and `select(term, estimate, conf.low, conf.high)`. Render. - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-3} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r tidymodels-3-test} -#| echo: true -linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area, data = ames) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) -``` - -### - -`fit()` is where the mechanism meets real data — this is step three of Structure → Method → Fitting. The `Gr_Liv_Area` coefficient is about 111, with an interval well clear of zero: each additional square foot is associated with about $111 more in sale price. - -### Exercise 4 - -Edit the formula to add `Year_Built` as a second predictor. Render. - -``` -linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-4} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r tidymodels-4-test} -#| echo: true -linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) -``` - -### - -Both intervals exclude zero — by this heuristic, `Year_Built` earns its place. Each additional year a house was built later is associated with about $1,087 more in sale price, holding living area fixed. - -### Exercise 5 - -What about overall condition? Edit the formula to add `Overall_Cond`. Render. - -``` -linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built + Overall_Cond, data = ames) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r tidymodels-5-test} -#| echo: true -linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built + Overall_Cond, data = ames) |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) -``` - -### - -`Overall_Cond` is categorical, so it produces one row per level — and the intervals are mixed: some levels clearly exclude zero, others don't. This is exactly the "many possible DGMs" problem: there's no automatic rule that resolves it. For this tutorial, we'll make the judgment call to keep things simple and go back to the two-variable model as our final DGM. - -### Exercise 6 - -Edit your working chunk one last time: remove `Overall_Cond`, remove the `tidy()`/`select()` steps, and assign the result to `mod_linear`. Print it. Render. - -``` -mod_linear <- linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames) -mod_linear -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r tidymodels-6} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r tidymodels-6-test} -#| echo: true -mod_linear <- linear_reg() |> - set_engine("lm") |> - fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames) -mod_linear -``` - -### - -From here on, you never touch this formula again. `mod_linear` is the fixed, saved model every later exercise addresses — the candidates you discarded never got a permanent name at all, which was deliberate: there are always more DGMs you *could* fit than you actually will, and we never know for certain which one is true. - -### Exercise 7 - -Add `#| cache: true` to the top of your `mod_linear` chunk, and add `*_cache/` to your `.gitignore`. In the R Terminal, run: - -``` -show_file(".gitignore") -``` - -CP/CR. - -```{r tidymodels-7} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 4) -``` - -### - -``` -analysis_files -*_cache/ -``` - -### - -Structure → Method → Fitting. You have now completed all three steps. `mod_linear` is your fitted DGM — the machine every later question will be addressed to. - -### Exercise 8 - -This exercise assumes `mod_linear` is exactly the model you assigned in Exercise 6. In a new code chunk, run the following to see `mod_linear`'s coefficients at a readable precision. Render. - -``` -mod_linear |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 0))) -``` - -As you can see, the coefficient on `Gr_Liv_Area` is about 96. What does that mean? - -```{r tidymodels-8} -question_text(NULL, - answer("For houses built in the same year, each additional square foot of above-ground living area is associated with an expected increase in sale price of $96, with a 95% confidence interval of $93 to $99.", correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -```{r tidymodels-8-test} -#| echo: true -mod_linear |> - tidy(conf.int = TRUE) |> - select(term, estimate, conf.low, conf.high) |> - mutate(across(where(is.numeric), \(x) round(x, 0))) -``` - -### - -Here's the answer: - -> For houses built in the same year, each additional square foot of above-ground living area is associated with an expected increase in sale price of $96, with a 95% confidence interval of $93 to $99. - -### - -This answer describes an expected, average change rather than what happens to any one specific house, and it holds `Year_Built` fixed while describing the change in `Gr_Liv_Area` — we never claim that adding square footage *causes* a $96 increase, only that houses differing by one square foot differ in expected price by about that much. - -### Exercise 9 - -Using the same table from Exercise 8: the coefficient on `Year_Built` is about 1,087. What does that mean? - -```{r tidymodels-9} -question_text(NULL, - answer("Holding above-ground living area constant, each additional year in a house's construction year is associated with an expected increase in sale price of $1,087, with a 95% confidence interval of $1,029 to $1,145.", correct = TRUE), - allow_retry = FALSE, - incorrect = NULL, - rows = 6) -``` - -### - -Here's the answer: - -> Holding above-ground living area constant, each additional year in a house's construction year is associated with an expected increase in sale price of $1,087, with a 95% confidence interval of $1,029 to $1,145. - -### - -Use that same pattern — expected value, other variables held fixed, no causal claim — whenever you're asked to interpret a coefficient in any model. The fitted DGM is not the answer to anything by itself; it is the machine we built so we can answer questions. - -### Exercise 10 - -Add the equation for `mod_linear` as a LaTeX math block above the model chunk: - -$$\widehat{\text{Sale\_Price}} = -2{,}106{,}459 + 96 \times \text{Gr\_Liv\_Area} + 1{,}087 \times \text{Year\_Built}$$ - -Render. In the R Terminal, run `show_file("analysis.qmd")`. CP/CR just the equation line and the `mod_linear` chunk right below it. - -```{r tidymodels-10} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 20) -``` - -### - -The relevant part of your output should look like this: - -
$$\widehat{\text{Sale\_Price}} = -2{,}106{,}459 + 96 \times \text{Gr\_Liv\_Area} + 1{,}087 \times \text{Year\_Built}$$
-
-#| cache: true
-mod_linear <- linear_reg() |>
-  set_engine("lm") |>
-  fit(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames)
-mod_linear
-
- -### - -This equation IS the fitted DGM, written out in full — any prediction it produces is just plugging numbers into it. A 1,500 sq ft house built in 1980 gets: $-2{,}106{,}459 + 96 \times 1500 + 1{,}087 \times 1980 \approx \$189{,}801$. - -### Exercise 11 - -Commit and push, with a message like "Fit final linear regression DGM for Ames housing." In the bash Terminal, run: - -``` -git log -1 -``` - -CP/CR. - -```{r tidymodels-11} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -``` -commit b2c3d4e5f6a7... -Author: Var Kurapati -Date: Mon Jul 27 2026 - Fit final linear regression DGM for Ames housing -``` - -### - -`mod_linear` is not the answer to anything by itself — it's the machine we built so we can answer questions, which is exactly what the next section does with it. - -## Marginal Effects -### - -We built the DGM with *Courage*. Now we practice *Temperance*: interpreting it, and using it — with humility — to answer the question we started with. - -Interpreting a model with **marginaleffects** almost always comes down to one distinction: **marginal** vs. **conditional**. So every time you see `by`, you are averaging over real houses, and every time you see `condition`, you are looking at a hypothetical. - -If Tidymodels was Structure → Method → Fitting, this section is **Question → Units → Answer**: pose the question as a quantity, specify which units you're asking about, then run the machine. - -In this section we cover five functions from [**marginaleffects**](https://marginaleffects.com/): `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()`. Keep asking yourself, for every function call in this section, which one you're looking at. - -### Exercise 1 - -In a new code chunk, run `predictions()` for one hypothetical house: 1,500 sq ft, built in 1980. Render. - -``` -predictions(mod_linear, - newdata = data.frame(Gr_Liv_Area = 1500, Year_Built = 1980)) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-1} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 6) -``` - -### - -```{r predictions-1-test} -#| echo: true -predictions(mod_linear, - newdata = data.frame(Gr_Liv_Area = 1500, Year_Built = 1980)) -``` - -### - -The fitted DGM predicts a price for this house — check the `estimate` column. It should be close to what the equation gives directly: $-2{,}106{,}459 + 96 \times 1500 + 1{,}087 \times 1980 \approx \$189{,}801$. This is a **conditional** prediction — nothing in Ames is actually this exact house. - -### Exercise 2 - -Update your chunk to compare two hypothetical houses that differ *only* in living area — 1,500 vs. 2,500 sq ft — with `Year_Built` held identical. Render. - -``` -predictions(mod_linear, - newdata = data.frame( - Gr_Liv_Area = c(1500, 2500), - Year_Built = c(1980, 1980))) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r predictions-2-test} -#| echo: true -predictions(mod_linear, - newdata = data.frame( - Gr_Liv_Area = c(1500, 2500), - Year_Built = c(1980, 1980))) -``` - -### - -The larger house is predicted to sell for about $96,000 more — one thousand extra square feet at $96 per square foot. That is what a coefficient of 96 *means* in practice. This is the counterfactual question: same age, different size, same machine asked twice. - -### Exercise 3 - -Update your chunk to predict the sale price for the first actual house in the dataset, and separately check its actual sale price. Render. - -``` -predictions(mod_linear, newdata = ames[1, ]) -ames$Sale_Price[1] -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-3} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r predictions-3-test} -#| echo: true -predictions(mod_linear, newdata = ames[1, ]) -ames$Sale_Price[1] -``` - -### - -The model's predicted price and the actual sale price won't match exactly. That gap — actual minus predicted — is the **residual**: the part of this house's price the model couldn't explain from just its size and age. No DGM is perfect; the residual is its honest accounting of what it missed. - -### Exercise 4 - -Update your chunk to call `predictions(mod_linear)` with no `newdata` at all, piped to `as_tibble()`, and print `nrow()` of the result. Render. - -``` -predictions(mod_linear) |> as_tibble() |> nrow() -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-4} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 6) -``` - -### - -```{r predictions-4-test} -#| echo: true -predictions(mod_linear) |> as_tibble() |> nrow() -``` - -### - -Without `newdata`, `predictions()` defaults to running the model on every row it was fit on — one prediction per real house. That number should match the 2,930 houses mentioned back in the Introduction. `as_tibble()` turns the result into something you can wrangle with ordinary **dplyr** verbs. - -### Exercise 5 - -Update your chunk to use `avg_predictions()` with `by = "Bldg_Type"` to get the average predicted sale price for each building type. Render. - -``` -avg_predictions(mod_linear, by = "Bldg_Type") -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -```{r predictions-5-test} -#| echo: true -avg_predictions(mod_linear, by = "Bldg_Type") -``` - -### - -`by` is always the same operation: run the model on every real house, then average within groups. This is a **marginal** prediction — it uses the actual distribution of houses in each building type, not a hypothetical one. - -### Exercise 6 - -Now build a **conditional** prediction using `datagrid()`: three hypothetical houses at 1,000, 2,000, and 3,000 sq ft, with `Year_Built` held at its typical value automatically. Render. - -``` -predictions(mod_linear, - newdata = datagrid(Gr_Liv_Area = c(1000, 2000, 3000))) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-6} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r predictions-6-test} -#| echo: true -predictions(mod_linear, - newdata = datagrid(Gr_Liv_Area = c(1000, 2000, 3000))) -``` - -### - -`datagrid()` fixes every variable you don't mention at a representative value — `Year_Built` lands at its mean — and varies only the ones you specify. This is the **conditional** counterpart to `avg_predictions()` with a `by` argument: no real house needs to exist at any of these three sizes. - -### Exercise 7 - -A prediction says what the model expects; a **comparison** says how much that expectation changes when a predictor changes. In a new code chunk, use `comparisons()` to compute the predicted change in sale price when `Gr_Liv_Area` increases by 100 sq ft, holding each house's other measurements at their own real values. Render. - -``` -cmp <- comparisons(mod_linear, - variables = list(Gr_Liv_Area = 100)) -nrow(cmp) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-7} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 6) -``` - -### - -```{r predictions-7-test} -#| echo: true -cmp <- comparisons(mod_linear, - variables = list(Gr_Liv_Area = 100)) -nrow(cmp) -``` - -### - -2,930 rows — one per real house. Each row is the predicted change in sale price for one specific house if its living area were 100 sq ft larger, everything else held fixed. For a linear model this number is the same for every house — but for logistic or ordered models it varies by unit, which is where `comparisons()` becomes essential. - -### Exercise 8 - -Update your chunk to use `avg_comparisons()` instead, to see the average predicted change across all real houses. Render. - -``` -avg_comparisons(mod_linear, - variables = list(Gr_Liv_Area = 100)) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-8} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r predictions-8-test} -#| echo: true -avg_comparisons(mod_linear, - variables = list(Gr_Liv_Area = 100)) -``` - -### - -About $9,600 per 100 sq ft — consistent with the $96 per sq ft coefficient from Tidymodels, translated into units a person can use. This is the coefficient from the model made concrete: not a log-odds or a slope, but a dollar amount for a real-world change. - -### Exercise 9 - -Update your chunk to ask the same question about age: the average predicted change in sale price if every house were built one year later. Render. - -``` -avg_comparisons(mod_linear, - variables = list(Year_Built = 1)) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-9} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 8) -``` - -### - -```{r predictions-9-test} -#| echo: true -avg_comparisons(mod_linear, - variables = list(Year_Built = 1)) -``` - -### - -About $1,087 per year — consistent with the coefficient from Tidymodels. `plot_comparisons()` will draw this difference visually: not two groups of different houses, but the same houses imagined in two worlds, one year apart. - -### Exercise 10 - -Ask your AI assistant to write code using `plot_predictions(mod_linear, condition = "Gr_Liv_Area", draw = FALSE)` piped into `ggplot()` to visualize how predicted sale price varies with `Gr_Liv_Area`, with a shaded confidence band, a title, axis labels, and a dollar-formatted y-axis. Paste the code into a new chunk. Render. - -``` -plot_predictions(mod_linear, condition = "Gr_Liv_Area", draw = FALSE) |> - ggplot(aes(x = Gr_Liv_Area, y = estimate)) + - geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) + - geom_line() + - labs( - title = "Predicted sale price increases with living area", - x = "Above-ground living area (sq ft)", - y = "Predicted sale price ($)" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-10} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 10) -``` - -### - -```{r predictions-10-test} -#| echo: true -plot_predictions(mod_linear, condition = "Gr_Liv_Area", draw = FALSE) |> - ggplot(aes(x = Gr_Liv_Area, y = estimate)) + - geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) + - geom_line() + - labs( - title = "Predicted sale price increases with living area", - x = "Above-ground living area (sq ft)", - y = "Predicted sale price ($)" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -### - -`draw = FALSE` returns the raw numbers behind the plot — `estimate`, `conf.low`, and `conf.high` — so you can build the chart yourself with full control over what it looks like. With `condition = "Gr_Liv_Area"`, the machine runs on hypothetical houses with `Year_Built` held at its mean. - -### Exercise 11 - -Ask your AI assistant to extend that plot to show predicted sale price across both `Gr_Liv_Area` and `Year_Built`, using `condition = c("Gr_Liv_Area", "Year_Built")` with `draw = FALSE`, piped into `ggplot()` — color the lines by `Year_Built`. Paste its code into a new chunk. Render. - -``` -plot_predictions(mod_linear, - condition = c("Gr_Liv_Area", "Year_Built"), draw = FALSE) |> - ggplot(aes(x = Gr_Liv_Area, y = estimate, color = factor(Year_Built))) + - geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(Year_Built)), - alpha = 0.15, color = NA) + - geom_line() + - labs( - title = "Predicted sale price by size and age", - x = "Above-ground living area (sq ft)", - y = "Predicted sale price ($)", - color = "Year Built", - fill = "Year Built" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-11} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 10) -``` - -### - -```{r predictions-11-test} -#| echo: true -plot_predictions(mod_linear, - condition = c("Gr_Liv_Area", "Year_Built"), draw = FALSE) |> - ggplot(aes(x = Gr_Liv_Area, y = estimate, color = factor(Year_Built))) + - geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(Year_Built)), - alpha = 0.15, color = NA) + - geom_line() + - labs( - title = "Predicted sale price by size and age", - x = "Above-ground living area (sq ft)", - y = "Predicted sale price ($)", - color = "Year Built", - fill = "Year Built" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -### - -When two variables are passed to `condition`, **marginaleffects** varies the first continuously and shows the second at representative values — that's what `factor(Year_Built)` groups into separate colored lines. Newer houses command higher prices at every size — the lines are parallel because the model has no interaction term. - -### Exercise 12 - -Now plot a **marginal** quantity instead: use `plot_predictions()` with `by = "Bldg_Type"` and `draw = FALSE`, piped into `ggplot()`, to visualize average predicted sale price per building type. Render. - -``` -plot_predictions(mod_linear, by = "Bldg_Type", draw = FALSE) |> - ggplot(aes(x = Bldg_Type, y = estimate)) + - geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + - labs( - title = "Average predicted sale price by building type", - x = NULL, - y = "Average predicted sale price ($)" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. - -```{r predictions-12} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 10) -``` - -### - -```{r predictions-12-test} -#| echo: true -plot_predictions(mod_linear, by = "Bldg_Type", draw = FALSE) |> - ggplot(aes(x = Bldg_Type, y = estimate)) + - geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + - labs( - title = "Average predicted sale price by building type", - x = NULL, - y = "Average predicted sale price ($)" - ) + - scale_y_continuous(labels = scales::dollar) -``` - -### - -Same function, same package — but passing `by` instead of `condition` switches `plot_predictions()` from showing a hypothetical relationship to showing an average over real houses, grouped by building type. `condition` is a conditional, counterfactual view; `by` is a marginal, descriptive one — the same numbers from the `avg_predictions(by = "Bldg_Type")` table earlier, drawn as a picture. - -### Exercise 13 - -Commit and push, with a message like "Add marginaleffects predictions, comparisons, and plots." In the bash Terminal, run: - -``` -git log -1 -``` - -CP/CR. - -```{r predictions-13} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 5) -``` - -### - -``` -commit c3d4e5f6a7b8... -Author: Var Kurapati -Date: Mon Jul 27 2026 - Add marginaleffects predictions, comparisons, and plots -``` - -### - -You have now used `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()` on a linear regression model, and drawn the marginal/conditional line clearly in both table and plot form. The same five functions work identically on logistic, multinomial, and ordered models — only the number of rows per case changes. - -## Summary -### - -You have built an analysis page that fits a linear regression model and visualizes its predictions. Now wire it into a navigable Quarto website and publish it. - -### Exercise 1 - -Create a `_quarto.yml` file at the top level of your repo with the following contents: - -``` -project: - type: website - -website: - title: "Ames Housing" - navbar: - left: - - href: analysis.qmd - text: Analysis - - href: about.qmd - text: About - -format: - html: - theme: cosmo -``` - -In the bash Terminal, run: - -``` -quarto render -ls -``` - -In the R Terminal, run `show_file("_quarto.yml")`. CP/CR. - -```{r summary-1} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 15) -``` - -### - -
project:
-  type: website
-
-website:
-  title: "Ames Housing"
-  navbar:
-    left:
-      - href: analysis.qmd
-        text: Analysis
-      - href: about.qmd
-        text: About
-
-format:
-  html:
-    theme: cosmo
-
- -### - -`_quarto.yml` turns a folder of `.qmd` files into a navigable website. Without it, each page renders independently. - -### Exercise 2 - -Create `about.qmd` with a YAML title of `"About"`, your name, and today's date. Also clean up `analysis.qmd` so it contains only: a setup chunk, the final `mod_linear` model chunk with its equation, and the final three `plot_predictions()` visualizations from the Marginal Effects section. Remove all intermediate output. Render. - -In the R Terminal, run `show_file("about.qmd")`, then `show_file("analysis.qmd")`. CP/CR both. - -```{r summary-2} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 12) -``` - -### - -``` ---- -title: "About" ---- - -Var Kurapati, July 2026. -``` - -### - -Your `analysis.qmd` should now have exactly five code chunks — setup, the `mod_linear` model, and the three `plot_predictions()` calls — plus one math block (the equation). - -### - -An About page tells visitors who made the site and when. - -### Exercise 3 - -Before publishing, preview the site locally. In the bash Terminal, run: - -``` -quarto preview -``` - -Click through the navbar in the browser preview and confirm both pages load with no errors. In the bash Terminal, stop the preview and CP/CR the first few lines of terminal output from the `quarto preview` command. - -```{r summary-3} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 6) -``` - -### - -``` -Preparing to preview -Watching files for changes -GET: / -GET: /analysis.html -GET: /about.html -Browse at http://localhost:3417/ -``` - -### - -`quarto preview` builds the site and opens a live-reloading local server — catching broken links or missing pages before they go live is much cheaper than catching them after publishing. - -### Exercise 4 - -Publish the website to GitHub Pages. In the bash Terminal, run: - -``` -quarto publish gh-pages -``` - -Copy/paste the resulting URL below. - -```{r summary-4} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 3) -``` - -### - -`quarto publish gh-pages` without a filename publishes the whole website. Both pages are now live, linked by the navigation bar. - -### Exercise 5 - -Commit and push any remaining changes. Copy/paste the URL to your GitHub repo. - -```{r summary-5} -question_text(NULL, - answer(NULL, correct = TRUE), - allow_retry = TRUE, - try_again_button = "Edit Answer", - incorrect = NULL, - rows = 3) -``` - -### - -The same pattern you used here — `linear_reg() |> set_engine() |> fit()` followed by `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()`, always keeping marginal and conditional straight — applies to logistic, multinomial, and ordered regression too. The interface never changes; only the structure does. - -```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} -``` diff --git a/inst/tutorials/dgm-logistic/tutorial.Rmd b/inst/tutorials/dgm-logistic/tutorial.Rmd index 4ffd590..926575a 100644 --- a/inst/tutorials/dgm-logistic/tutorial.Rmd +++ b/inst/tutorials/dgm-logistic/tutorial.Rmd @@ -610,7 +610,7 @@ Interpreting a model with **marginaleffects** almost always comes down to one di If Tidymodels was Structure → Method → Fitting, this section is **Question → Units → Answer**: pose the question as a quantity, specify which units you're asking about, then run the machine. -In this section we cover five functions from [**marginaleffects**](https://marginaleffects.com/): `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()`. The one thing that changes for a logistic model: every prediction comes back as one row per outcome category, so you'll filter to `group == "Yes"` to keep just the probability of leaving. +In this section we cover four functions from [**marginaleffects**](https://marginaleffects.com/): `predictions()`, `avg_predictions()`, `comparisons()`, and `avg_comparisons()` --- and build every plot directly from their output with `ggplot()`. The one thing that changes for a logistic model: a prediction comes back as one row per outcome category, so the printed tables filter to `group == "Yes"` to keep just the probability of leaving. The plotting exercises instead run these functions on `extract_fit_engine(mod_logistic)` --- the underlying `glm` --- whose `type = "response"` scale is that probability directly, with no `group` column to filter. ### Exercise 1 @@ -768,11 +768,10 @@ This employee's outcome is already known — so why predict it? Because the fitt ### Exercise 5 -Update your chunk to use `avg_predictions()` with `by = "OverTime"` and `type = "prob"`, filtered to `group == "Yes"`. Render. +Update your chunk to use `avg_predictions()` with `by = "OverTime"`, run on `extract_fit_engine(mod_logistic)` --- the underlying `glm`, whose predictions are the probability of leaving directly, so there is no `group` column to filter. Render. ``` -avg_predictions(mod_logistic, by = "OverTime", type = "prob") |> - filter(group == "Yes") +avg_predictions(extract_fit_engine(mod_logistic), by = "OverTime") ``` In the R Terminal, run `show_file("analysis.qmd", chunk = "Last")`. CP/CR. @@ -790,13 +789,12 @@ question_text(NULL, ```{r predictions-5-test} #| echo: true -avg_predictions(mod_logistic, by = "OverTime", type = "prob") |> - filter(group == "Yes") +avg_predictions(extract_fit_engine(mod_logistic), by = "OverTime") ``` ### -This is a **marginal** prediction — run the model on every real employee, then average within overtime groups. The gap between the two rows is the overtime effect, averaged over the actual income and tenure distribution of real employees in each group. +This is a **marginal** prediction — run the model on every real employee, then average within overtime groups. The gap between the two rows is the overtime effect, averaged over the actual income and tenure distribution of real employees in each group. We call `avg_predictions()` on `extract_fit_engine(mod_logistic)`, the raw `glm` behind the tidymodels wrapper: **marginaleffects** returns the grouping column and the response-scale probability more cleanly from the engine object, which matters once we start plotting. ### Exercise 6 @@ -944,13 +942,17 @@ About −1.4 percentage points per $1,000 — a modest effect compared to the 15 ### Exercise 10 -Ask your AI assistant to write `plot_predictions()` code that visualizes how predicted attrition probability varies with `MonthlyIncome`, using `mod_logistic`, `type = "prob"`, and `condition = "MonthlyIncome"` with `draw = FALSE`, piped into a `filter(group == "Yes")` step then `ggplot()`. Render. +Ask your AI assistant to write code that visualizes how predicted attrition probability varies with `MonthlyIncome`. Build a grid of hypothetical employees --- `MonthlyIncome` swept across its range, `YearsAtCompany` at its mean, `OverTime` set to `"No"` --- pass it to `predictions()` on `extract_fit_engine(mod_logistic)` with `type = "response"`, then `ggplot()`. Render. ``` -plot_predictions(mod_logistic, - condition = "MonthlyIncome", - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +grid <- tibble( + MonthlyIncome = seq(min(attrition$MonthlyIncome), + max(attrition$MonthlyIncome), length.out = 100), + YearsAtCompany = mean(attrition$YearsAtCompany), + OverTime = "No" +) + +predictions(extract_fit_engine(mod_logistic), newdata = grid, type = "response") |> ggplot(aes(x = MonthlyIncome, y = estimate)) + geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) + geom_line() + @@ -978,10 +980,14 @@ question_text(NULL, ```{r predictions-10-test} #| echo: true -plot_predictions(mod_logistic, - condition = "MonthlyIncome", - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +grid <- tibble( + MonthlyIncome = seq(min(attrition$MonthlyIncome), + max(attrition$MonthlyIncome), length.out = 100), + YearsAtCompany = mean(attrition$YearsAtCompany), + OverTime = "No" +) + +predictions(extract_fit_engine(mod_logistic), newdata = grid, type = "response") |> ggplot(aes(x = MonthlyIncome, y = estimate)) + geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) + geom_line() + @@ -996,17 +1002,21 @@ plot_predictions(mod_logistic, ### -`plot_predictions()` draws the machine's answers. With `condition = "MonthlyIncome"`, the machine runs on hypothetical employees with `YearsAtCompany` and `OverTime` held at their typical values. The curve bends rather than running in a straight line — probabilities are squeezed between 0 and 1, so the same change in log-odds produces a smaller change in probability at the extremes. +The grid holds `YearsAtCompany` and `OverTime` fixed while `MonthlyIncome` sweeps its range; `predictions()` on `extract_fit_engine(mod_logistic)` --- the raw `glm` --- runs the machine on each hypothetical employee and returns the probability of leaving directly on the response scale. The curve bends rather than running in a straight line — probabilities are squeezed between 0 and 1, so the same change in log-odds produces a smaller change in probability at the extremes. ### Exercise 11 -Ask your AI assistant to extend that plot to show predicted attrition probability across both `MonthlyIncome` and `OverTime`, using `condition = c("MonthlyIncome", "OverTime")` — color the lines by `OverTime`. Render. +Ask your AI assistant to extend that plot to show predicted attrition probability across both `MonthlyIncome` and `OverTime`. Rebuild the grid with `expand_grid()` so the income sweep is crossed with both `OverTime` values (keep `YearsAtCompany` at its mean), and color the lines by `OverTime`. Render. ``` -plot_predictions(mod_logistic, - condition = c("MonthlyIncome", "OverTime"), - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +grid <- expand_grid( + MonthlyIncome = seq(min(attrition$MonthlyIncome), + max(attrition$MonthlyIncome), length.out = 100), + OverTime = c("No", "Yes") +) |> + mutate(YearsAtCompany = mean(attrition$YearsAtCompany)) + +predictions(extract_fit_engine(mod_logistic), newdata = grid, type = "response") |> ggplot(aes(x = MonthlyIncome, y = estimate, color = OverTime)) + geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = OverTime), alpha = 0.15, color = NA) + @@ -1035,10 +1045,14 @@ question_text(NULL, ```{r predictions-11-test} #| echo: true -plot_predictions(mod_logistic, - condition = c("MonthlyIncome", "OverTime"), - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +grid <- expand_grid( + MonthlyIncome = seq(min(attrition$MonthlyIncome), + max(attrition$MonthlyIncome), length.out = 100), + OverTime = c("No", "Yes") +) |> + mutate(YearsAtCompany = mean(attrition$YearsAtCompany)) + +predictions(extract_fit_engine(mod_logistic), newdata = grid, type = "response") |> ggplot(aes(x = MonthlyIncome, y = estimate, color = OverTime)) + geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = OverTime), alpha = 0.15, color = NA) + @@ -1054,17 +1068,14 @@ plot_predictions(mod_logistic, ### -When two variables are passed to `condition`, **marginaleffects** varies the first continuously and shows the second at its distinct values — that's what `color = OverTime` splits into two separate curves. The gap between the curves is the overtime effect, which the `avg_comparisons()` exercises measured as about 15 percentage points. +`expand_grid()` crosses the income sweep with both overtime values, so `predictions()` returns one curve per overtime status — that's what `color = OverTime` splits apart. The gap between the curves is the overtime effect, which the `avg_comparisons()` exercises measured as about 15 percentage points. ### Exercise 12 -Now plot a **marginal** quantity instead: use `plot_predictions()` with `by = "OverTime"`, `type = "prob"`, and `draw = FALSE`, filtered to `group == "Yes"`, to visualize average predicted attrition probability by overtime group. Render. +Now plot a **marginal** quantity instead: use `avg_predictions()` with `by = "OverTime"` on `extract_fit_engine(mod_logistic)`, as in Exercise 5, piped into `ggplot()`, to visualize average predicted attrition probability by overtime group. Render. ``` -plot_predictions(mod_logistic, - by = "OverTime", - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +avg_predictions(extract_fit_engine(mod_logistic), by = "OverTime") |> ggplot(aes(x = OverTime, y = estimate)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + labs( @@ -1090,10 +1101,7 @@ question_text(NULL, ```{r predictions-12-test} #| echo: true -plot_predictions(mod_logistic, - by = "OverTime", - type = "prob", draw = FALSE) |> - filter(group == "Yes") |> +avg_predictions(extract_fit_engine(mod_logistic), by = "OverTime") |> ggplot(aes(x = OverTime, y = estimate)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + labs( @@ -1106,7 +1114,7 @@ plot_predictions(mod_logistic, ### -Same function, same package — but passing `by` instead of `condition` switches from a hypothetical relationship to an average over real employees, grouped by overtime status. `condition` is a conditional, counterfactual view; `by` is a marginal, descriptive one — the same two numbers `avg_predictions()` gave you earlier, drawn as a picture. +`avg_predictions()` with `by` runs the model on every real employee and averages within groups — a marginal, descriptive view, in contrast to the conditional, counterfactual grids of the previous two exercises. These are the same two numbers you computed in Exercise 5, now drawn as a picture. ### Exercise 13 @@ -1138,7 +1146,7 @@ Date: Mon Aug 03 2026 ### -You have now used `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()` on a logistic regression model. The same five functions work identically on multinomial and ordered models — except there, filtering by `group` becomes more interesting, since there are more than two categories to choose from. +You have now used `predictions()`, `avg_predictions()`, `comparisons()`, and `avg_comparisons()` on a logistic regression model, building each plot from their output with `ggplot()`. The same functions work identically on multinomial and ordered models — except there, filtering by `group` becomes more interesting, since there are more than two categories to choose from. ## Summary ### @@ -1210,7 +1218,7 @@ format: ### Exercise 2 -Create `about.qmd` with a YAML title of `"About"`, your name, and today's date. Clean up `analysis.qmd` so it contains only: a setup chunk, the final `mod_logistic` model chunk with its equation, and the final three `plot_predictions()` visualizations. Remove all intermediate output. Render. +Create `about.qmd` with a YAML title of `"About"`, your name, and today's date. Clean up `analysis.qmd` so it contains only: a setup chunk, the final `mod_logistic` model chunk with its equation, and the final three prediction visualizations (the two `predictions()` plots and the `avg_predictions()` plot). Remove all intermediate output. Render. In the R Terminal, run `show_file("about.qmd")`, then `show_file("analysis.qmd")`. CP/CR both. @@ -1235,7 +1243,7 @@ Var Kurapati, August 2026. ### -Your `analysis.qmd` should now have exactly five code chunks — setup, the `mod_logistic` model, and the three `plot_predictions()` calls — plus one math block. +Your `analysis.qmd` should now have exactly five code chunks — setup, the `mod_logistic` model, and the three prediction plots — plus one math block. ### @@ -1313,7 +1321,7 @@ question_text(NULL, ### -The same pattern you used here — `logistic_reg() |> set_engine() |> fit()` followed by `predictions()`, `avg_predictions()`, `comparisons()`, `avg_comparisons()`, and `plot_predictions()`, always keeping marginal and conditional straight — applies to multinomial and ordered regression too. Structure → Method → Fitting built the DGM. Question → Units → Answer ran it. Every tutorial in this series is just those six words, pointed at a different kind of outcome. +The same pattern you used here — `logistic_reg() |> set_engine() |> fit()` followed by `predictions()`, `avg_predictions()`, `comparisons()`, and `avg_comparisons()`, always keeping marginal and conditional straight — applies to multinomial and ordered regression too. Structure → Method → Fitting built the DGM. Question → Units → Answer ran it. Every tutorial in this series is just those six words, pointed at a different kind of outcome. ```{r download-answers, child = system.file("child_documents/download_answers.Rmd", package = "tutorial.helpers")} ```