Probability Practice with NSDUH · Cross-Tabs and Distributions

Lab · Module 6 · Wed Sep 23

Welcome to the Module 6 lab

The M06 module and pre-study cover both probability rules and distributions, and today’s lab gives you practice with both:

  • The probability-rules half gave you the rules of probability — marginals, joints, conditionals, the addition rule. Part 1 of today’s lab applies those rules to a real clinical question using a 2×2 cross-tabulation.
  • The distributions half introduced the normal distribution and the pnorm() / qnorm() machinery. Part 2 adds the discipline those examples assumed: estimate the model from a sample, check that it fits, then use it.

Two kinds of box — and only one goes in your notebook

What it looks like What to do
Check yourself a box to type in or options to click, with a Check your thinking note underneath Answer it right here on the page. It marks itself. Nothing goes in your notebook.
Checkpoint a heather box at the end of each step The one box that matters for your grade. It names the section your notebook should now hold, gives you the number to check your work against, and — three times in this lab — asks for a sentence or two of prose.

Everything the checkpoints ask you to write is listed again in the Double check box at the end.

The substantive thread runs through both parts:

Among adolescent females in the 2019 NSDUH analytic sample, how often do a past-year major depressive episode (MDE) and a past-year substance use disorder (SUD) co-occur — and how does conditioning on one change the observed probability of the other? Then, across all MDE-positive adolescents, how does depression’s interference distribute?

You’ll work with 2019 data from the National Survey on Drug Use and Health (NSDUH) — the same file you met in the M05 Module and pre-study. Because this lab uses raw, unweighted counts, we’ll treat every number as a descriptive sample probability — a summary of these respondents, not a survey-weighted estimate for all U.S. adolescents (those require additional tools we won’t use today).

How the time works. You’ll do the analysis — Part 1 and Part 2 — in class, then we’ll pause for a short group debrief. The write-up (the written interpretations each step asks for, the final render, and submitting to Canvas) is yours to finish at home. Draft those sentences roughly in the room while the output is in front of you, then tighten them at home — don’t rush the prose in class, it’s the part that takes the most thought.

What you’ll leave with

A rendered m06_lab.html report in programs/ containing:

  • The MDE × SUD cross-tab in four forms — counts, then cell, row, and column percentages — each read as a specific probability (Part 1).
  • A composite severity score built row-wise from four domain items (Part 2).
  • A fitted normal model for that score, estimated from the sample.
  • Model-vs-empirical probability statements from pnorm() and qnorm(), compared against what the data actually show.
  • Prose interpretation at three points, listed in the Double check box at the end.

nsduh_2019 · 13,397 observations · NSDUH 2019 wave · data/nsduh_2019.Rds

The 2019 wave of the survey you met in the M05 Module and pre-study. Part 1 uses three categorical variables; Part 2 uses four severity items that get averaged into a single composite score.

Part 1 variables (cross-tabs):

  • sex factor — Respondent’s sex as recorded by NSDUH
  • mde_pastyear factor — Whether the respondent experienced a major depressive episode in the past 12 months
  • substance_disorder factor — Whether the respondent met criteria for a past-year alcohol or illicit-drug use disorder (abuse or dependence)

Part 2 variables (distribution fitting):

  • severity_chores numeric — Sheehan Disability Scale rating of how much past-year depression interfered with household chores, where 0 is no interference and 10 is very severe interference
  • severity_work numeric — Sheehan Disability Scale rating of interference with school or work, 0 to 10
  • severity_family numeric — Sheehan Disability Scale rating of interference with family relationships, 0 to 10
  • severity_social numeric — Sheehan Disability Scale rating of interference with social life, 0 to 10

Full codebook for nsduh_2019 — worth opening once, because it explains why the severity items are missing for most rows and what the 0 on that scale means.

The dataset is already loaded into every sandbox chunk on this page under the name nsduh_2019, and Part 2’s derived objects (nsduh_mde_2019, mu_hat, sd_hat) are pre-built behind the scenes so the Part-2 chunks run independently of the Part-1 chunks. You’ll still build everything yourself in your RStudio notebook.

How this page works

Your real work happens in your RStudio notebook — the rendered HTML report is the deliverable, and the Steps below walk you through the sections you add to it. The sandbox chunks on this page run R live in the browser, using the same three-tab panel as the pre-studies:

  • ✍️ Your Code — an empty box you write in yourself
  • 💡 Hint — a nudge if you’re stuck
  • 👀 Spoiler — the full working code

Something has changed since Part 1. Through M05 the ✍️ tab handed you a partly written pipe with blanks to fill. From this lab on it is empty — except where a line is explicitly marked GIVEN, as in Step 8 — and a Targets list above each box tells you what to build — which object to create, which functions to reach for, and what to print to check yourself. Writing it from a blank box is harder, and that difficulty is the point: recalling a pipe is what makes it yours, and in M10–M12 there is no box at all — just your own notebook.

Work through the ✍️ tab first on every sandbox chunk. If you’ve been stuck for more than a few minutes, open 💡, and then 👀 if you need it — and keep moving. Getting through the whole lab matters more than getting any one chunk unaided.

Before you start

Two short setup steps before the two main parts of the lab.

Step 0 · Get set up

  1. Start in GitHub Desktop, before RStudio. Select PSY652_project, click Fetch origin, and click Pull origin if it appears. That is the pull half of the pull → edit → commit → push loop — the same loop you’ll close at the end of today’s lab.

  2. Now open RStudio by double-clicking PSY652_project.Rproj.

  3. In the Files pane, open programs/lab_template.qmd.

  4. Use File → Save As… and save a copy of the template inside programs/ with the name:

    m06_lab.qmd

  5. In the YAML header of your new file, update the title:

    title: "M06 Lab: NSDUH Cross-Tabs and Distributions"

  6. Delete the template’s three stock sections — # Setup, # Import Data (which reads a dataset called nhanes purely as a worked example), and # Glimpse. You are replacing them with your own structure in the next item.

  7. Sketch the report’s structure. Lay out all nine sections now and work down the list as the lab proceeds — you fill # Setup in a moment and # Import data in Step 1; every step after that fills one of the rest:

    # Setup
# Import data
# Cross-tabulation: MDE and SUD in females
# Probabilities from the cross-tab
# Subset to 2019 MDE-positive adolescents
# Build a composite severity score
# Describe the severity distribution
# Fit a normal model to severity_scale
# Probability statements from the fitted normal

    Each section should end up holding its code and a line or two of prose, so a reader can follow what you did without reading the code.

  8. If any line of that header or those chunks is a mystery, open lab_template_annotated.qmd — it sits beside the template in programs/ and is the same skeleton with a note on every YAML line, chunk option, and section. You will not need it every week, but it is the reference when something in the header stops making sense. From M07 on, the labs will point you back to this Step 0 rather than repeat it.

Insert your own R chunk before you carry any code over

The code boxes on this page run in the browser. Your own notebook runs on your copy of R, so anything you bring across has to sit inside an R code chunk that you create in your .qmd.

Insert the chunk first, the way you have since M02: click the green +C button at the top of the Source pane and choose R, or press Cmd+Option+I (Mac) / Ctrl+Alt+I (Windows).

Then give the chunk a label and put the code inside it:

#| label: glimpse-data

my_data |> glimpse()

Labels were optional in M03; from here on they are not. Now that you build the notebook yourself, put a short name on every chunk — glimpse-data, age-histogram. Quarto then names the chunk that failed instead of counting them, and the Outline button maps a long notebook for you.

If code lands outside a chunk, nothing runs. Quarto treats it as ordinary text and prints it in your report unexecuted. You usually find out several chunks later, when you get object 'xxx' not found naming a variable that looks perfectly fine — scroll up and check that the code creating it is genuinely inside a chunk.

Load the packages

Under your # Setup heading, add a chunk that loads the four packages this lab needs. Keep the label the template used:

#| label: setup

library(tidyverse)
library(here)
library(skimr)
library(gtsummary)

Then click Render — if m06_lab.html appears in programs/ with no red package errors, you are set up.

Step 1 · Import the data

In your notebook — this step fills your # Import data section.

Add a chunk labelled import-data that reads this week’s file:

nsduh_2019 <- read_rds(here("data", "nsduh_2019.Rds"))

Add a glimpse() line underneath it so you can see what came in, then run the chunk. A new object called nsduh_2019 should appear in your Environment tab.

Checkpoint 1 · Data imported

In your notebook — your # Import data section holds the import-data chunk and a glimpse() of what it read.

Check it worked — Your Environment pane shows an object called nsduh_2019 with 13,397 rows. If it isn’t there, confirm you named the object nsduh_2019 and pointed the path at data/nsduh_2019.Rds.

Part 1 · Discrete probability with cross-tabs

The substantive question for Part 1:

Among adolescent females in the 2019 NSDUH sample, how often does a past-year major depressive episode co-occur with a past-year substance use disorder? And does knowing someone has one update your belief that they also have the other?

Step 2 · A first cross-tabulation

In your notebook — this step fills your # Cross-tabulation: MDE and SUD in females section.

Check the levels before you filter them. filter() matches a factor level exactly, and a capital letter in the wrong place returns zero rows with no error and no warning. So ask the data rather than guessing — count() shows every level as R stores it:

Run it, then change mde_pastyear to sex and to substance_disorder to see those too. Get in the habit: count() a factor before you filter it. It takes three seconds and it is how you catch a level you misremembered, an unexpected NA, or a category with far fewer rows than you assumed.

Targets. The pipe you need does a few things in sequence — a common wrangling-then-table pattern you’ll reuse all semester:

  1. Start from nsduh_2019
  2. filter() to females (we’ll compare males later)
  3. select() only the two variables we need — that tightens the output and avoids accidental extra columns
  4. drop_na() to remove rows missing on either selected variable
  5. Pipe into tbl_cross() with mde_pastyear as rows and substance_disorder as columns

The five-step pattern above: filter() to sex == "Female" (capitalized level, in quotes), select() the two variables mde_pastyear and substance_disorder, drop_na(), then tbl_cross() with row = mde_pastyear and col = substance_disorder.

The table shows raw counts: how many females fell into each combination of MDE × SUD status. Read it carefully — every probability in the next step lives inside these four cells plus the row and column totals.

Checkpoint 2 · Cross-tab prints

In your notebook — your # Cross-tabulation: MDE and SUD in females section now holds this step’s code and output.

Check it worked — Your female-only cross-tab renders as a 2×2 table of counts, and the four interior cells sum to your female complete-case analytic n — the grand total in the table’s corner. Every probability in the next step is hiding somewhere in those numbers.

Step 3 · Read the probabilities off the table

In your notebook — this step fills your # Probabilities from the cross-tab section.

Every number the probability rules named — marginal, joint, conditional, union — is already in the 2×2 table you just built. You do not have to compute them by hand: tbl_cross() will do the arithmetic if you tell it which denominator you want.

The percent argument takes three values, and each one divides by a different total — which is what makes it answer a different question:

percent = Divides each cell by Gives you
"cell" the grand total joint probabilities inside the table, marginal probabilities at the margins
"row" that cell’s row total conditional probabilities given the row — P(SUD | MDE)
"column" that cell’s column total conditional probabilities given the column — P(MDE | SUD)

Targets. Rebuild the Step 2 pipe — start from nsduh_2019, filter() to females, select() the two variables, drop_na() — then add percent to the tbl_cross() call. Run it three times, changing only that one word, and keep all three tables in your notebook.

It is the Step 2 pipe with one argument added: percent = "cell", then "row", then "column". Nothing else changes.

With those three tables in front of you, every probability below is a number you can point to.

All three appeared in the M06 pre-study. Read each one aloud the way it is written in plain English:

Symbol Read aloud as Plain English Example
∩ “and” (intersection) both events happen to the same person P(MDE ∩ SUD) = “probability of MDE and SUD”
∪ “or” (union) at least one of the events happens P(MDE ∪ SUD) = “probability of MDE or SUD (or both)”
| “given” (conditional) restricted to a subset, then asks about the event P(SUD | MDE) = “probability of SUD given MDE”

Every probability statement you’ll read in a stats paper is some combination of these three.

Work through each of these

(a) Marginal probability of MDE: P(MDE).

On a cross-tab, marginals live in the row and column totals — the margins of the table.

\[P(\text{MDE}) = \frac{\text{females in the complete-case table who have MDE}}{\text{all females in the complete-case table}}\]

(That is the denominator for every formula in this step. Below we shorten it to “total females” for readability — but it always means the complete-case total.)

Talk it through

24% of adolescent females in this table had a past-year major depressive episode. Does that number surprise you? Thirty seconds with your neighbour — and say which direction it surprised you, because that is the more interesting half of the answer.

(b) Marginal probability of SUD: P(SUD).

Same idea as (a), now for substance use disorder. Pull it from the column totals.

Check yourself. Which condition is more common in this sample?

Roughly how many times more common is it? (one decimal place)

(c) Joint probability: P(MDE ∩ SUD), read “P of MDE and SUD”.

Joints live inside the table — in the four interior cells, not the margins.

\[P(\text{MDE} \cap \text{SUD}) = \frac{\text{total females with both}}{\text{total females}}\]

Check yourself. The joint can never exceed either marginal. Why must that be true?

Being in the overlap means having both — so every person counted in the joint is also counted in each marginal. A part cannot be bigger than the whole it sits inside.

That also tells you when the joint would equal a marginal: only if everyone with the rarer condition also had the other. Here, 50% of those with SUD also had MDE — high, but not everyone, so the joint stays below both.

(d) Union of two events — P(MDE ∪ SUD), read “P of MDE or SUD”.

The union has no single home on the table — you assemble it from three numbers you already have. Adding P(MDE) + P(SUD) would double-count the people who have both, so the addition rule subtracts them back out:

\[P(\text{MDE} \cup \text{SUD}) = P(\text{MDE}) + P(\text{SUD}) - P(\text{MDE} \cap \text{SUD})\]

Check yourself. Your union should come out only a little above P(MDE) on its own. Why so little?

The union adds only the people with SUD who don’t already have MDE. SUD is uncommon to begin with, and roughly half of those cases already sit inside the MDE group — so the union lands only 2.7 points above P(MDE). Overlap is what keeps a union small.

(e) Conditional probabilities — P(SUD | MDE) and P(MDE | SUD), read “P of SUD given MDE” and “P of MDE given SUD”.

A conditional restricts you to one row or one column and asks about the other variable within it — so the denominator is that row or column total, never the grand total.

Formal notation:

\[P(\text{SUD} \mid \text{MDE}) = \frac{P(\text{MDE} \cap \text{SUD})}{P(\text{MDE})} = \frac{\text{count with both}}{\text{total with MDE}}\]

\[P(\text{MDE} \mid \text{SUD}) = \frac{P(\text{MDE} \cap \text{SUD})}{P(\text{SUD})} = \frac{\text{count with both}}{\text{total with SUD}}\]

Same thing in plain English:

\[P(\text{SUD} \;\text{given}\; \text{MDE}) = \frac{P(\text{MDE} \;\text{and}\; \text{SUD})}{P(\text{MDE})}\]

\[P(\text{MDE} \;\text{given}\; \text{SUD}) = \frac{P(\text{MDE} \;\text{and}\; \text{SUD})}{P(\text{SUD})}\]

Notice that the joint probability (the numerator) is the same in both formulas — it’s the denominator that changes, because you’re restricting to a different subset each time.

Check yourself. Which of your two conditionals is larger?

The two are strikingly different — about 11% against roughly 50%. Knowing someone has a substance use disorder tells you far more about their depression status than the reverse.

And the reason is base rates, not a stronger link in one direction. Put the two fractions side by side and the only thing that moves is the denominator:

\[P(\text{SUD} \mid \text{MDE}) = \frac{170}{1,487} = 0.11 \qquad\text{versus}\qquad P(\text{MDE} \mid \text{SUD}) = \frac{170}{341} = 0.50\]

The numerator is the same 170 females in both — the overlap does not change depending on which way you read it. What changes is what you divide by: there are 1,487 females with MDE in this table but only 341 with an SUD. Dividing that shared overlap by the smaller SUD group gives the larger probability; dividing it by the bigger MDE group gives the smaller one. Neither number says the “link” runs more strongly one way — they are the same association, viewed from two differently sized starting groups.

This is exactly the asymmetry M06’s diagnostic-testing bonus turns into positive predictive value — and it is the single most common place people misread a conditional probability.

Checkpoint 3 · Six probabilities written

In your notebook — your # Probabilities from the cross-tab section now holds this step’s code and output.

Check it worked — Your notebook shows the cross-tab in all four forms — counts, cell %, row %, and column % — and your # Probabilities from the cross-tab section names all six probabilities in prose, each with the denominator it came from.

Write it up — the six probabilities as narrative text (not a code chunk), one line each, naming the denominator you divided by — then a short paragraph (3–5 sentences) tying them together. Some prompts to help:

  • Does conditioning on MDE meaningfully change the observed probability of SUD? Compare \(P(\text{SUD} \mid \text{MDE})\) with \(P(\text{SUD})\) — then do the reverse, comparing \(P(\text{MDE} \mid \text{SUD})\) with \(P(\text{MDE})\).
  • Which of the six numbers would you lead with if you were writing a one-paragraph abstract about MDE–SUD co-occurrence among adolescent females?

Part 2 · Continuous probability with the normal model

A question the pre-study didn’t ask

The M06 pre-study used the normal distribution to answer real probability questions — pnorm() for tail areas, qnorm() for cutoffs. Every example came with an assurance up front that the normal model was a reasonable fit for the variable at hand — and in the birth-weight Final Practice the parameters were simply handed to you as part of a stated model (\(\mu = 3254\) g, \(\sigma = 588\) g).

Real research often isn’t like that. You don’t know the true parameters — you estimate them from your sample and call the estimates \(\hat{\mu}\) and \(\hat{\sigma}\) (“mu-hat” and “sigma-hat”). You often don’t know the shape of the distribution either — so a probability read off a fitted model can be badly wrong, especially in the tails, which is usually where the clinically interesting decisions happen. Part 2 ends by checking the model’s answer against the data itself; the optional fit-checking bonus adds the two visual tools that catch a bad fit before you use it.

The analytic sample changes here — read this before Step 4

Part 1 was adolescent females only. Every probability you computed came from a table you had filtered with sex == "Female", because the question was about co-occurrence within one group.

Part 2 uses adolescents of all sexes. The only filter from here on is mde_pastyear == "Positive" — there is no sex filter in Part 2, and you should not add one. The question has changed: instead of asking how two conditions co-occur within a group, we are asking how depression severity is distributed among everyone who screened positive, and we want the largest sample we can get for estimating \(\hat{\mu}\) and \(\hat{\sigma}\).

So the row count you print in Step 4 will not be a subset of the female table from Part 1 — it is a differently defined group. If your Step 4 nrow() looks smaller than you expected, check that you have not carried the sex == "Female" filter forward out of habit.

Step 4 · Subset to MDE-positive adolescents

In your notebook — this step fills your # Subset to 2019 MDE-positive adolescents section.

The severity items are only measured for respondents who screened positive for past-year MDE. So Part 2 works on a stricter subset of nsduh_2019 than Part 1 did.

Use filter() to keep only the rows where mde_pastyear is "Positive", and select() the four severity items (plus year and mde_pastyear for record-keeping).

Targets.

  1. Start from nsduh_2019 and filter() to the rows where mde_pastyear is "Positive".
  2. select() year, mde_pastyear, and the four severity items.
  3. Assign the result to nsduh_mde_2019 — every later step uses that name, so match it exactly.
  4. Print nrow() and glimpse() on it to confirm the subset is the size and shape you expect before you build anything on top of it.

filter() to mde_pastyear == "Positive", then select() year, mde_pastyear, and starts_with("severity_") — the starts_with() helper grabs all four severity items without naming them one by one.

nsduh_mde_2019 should have 2,098 rows — that’s the count of 2019 adolescents who screened positive for past-year MDE.

Checkpoint 4 · MDE-positive subset built

In your notebook — your # Subset to 2019 MDE-positive adolescents section now holds this step’s code and output.

Check it worked — Your Environment pane shows nsduh_mde_2019 with 2,098 rows and the four severity_ columns kept. This stricter subset is the one every Part-2 calculation uses.

Step 5 · Build a composite severity score

In your notebook — this step fills your # Build a composite severity score section.

NSDUH measures depression’s interference across four life domains, one item each:

  • severity_chores — interference with home chores
  • severity_work — interference with school/work
  • severity_family — interference with family relationships
  • severity_social — interference with social life

Each item is rated on a 0–10 scale.

NSDUH already has a severity rule — and it is not the one we are about to build

You verified this yourself in M05. NSDUH scores role impairment from the maximum across the four domains, treating a rating of 7 or higher in any single domain as severe impairment. The binary result already sits in your data as mde_pastyear_severe, and in the M05 Module you rebuilt that rule by hand and confirmed it reproduces the column exactly.

What you are about to build is deliberately different. severity_scale is the mean across the four domains — a teaching index for practising row-wise scale construction, not NSDUH’s severity measure. The two rules answer different questions: a mean asks “how impaired was this adolescent across life domains, on average?”, while NSDUH’s maximum asks “was this adolescent severely impaired in at least one domain?” An adolescent reporting 10 on school/work and 0 elsewhere averages 2.5 — mild by our composite — but is severe by NSDUH’s.

Why a mean, for this lab? Because we need a variable with enough distinct values to make continuous-model fitting meaningful. The maximum of four integers is still an integer, so it takes only 11 possible values (0 through 10) — a chunky histogram no smooth curve fits usefully. The mean of four integers lands on quarter-points, giving 41 possible values, which is enough granularity to fit a normal model and see honestly where it succeeds and fails.

The M05 Module introduced the standard R pattern — rowMeans() combined with pick(). pick() selects a set of columns to operate on; rowMeans() averages across them. Together they compute the per-respondent mean in one line.

Targets.

  1. mutate() a new column severity_scale on nsduh_mde_2019 that is the row-wise mean of the four severity items.
  2. drop_na() on severity_scale, so a respondent missing any domain item is dropped rather than averaged from a partial set.
  3. Assign the result back to nsduh_mde_2019.
  4. Print the four domain items alongside the new composite for the first few rows, so you can check by eye that the average is right.

Build the composite with rowMeans(pick(starts_with("severity_"))) — pick() selects the four items, rowMeans() averages across them per row. The trailing drop_na() removes respondents missing an item (rowMeans() returns NA when any input is missing).

Scan the first few rows. Each severity_scale value is the average of the four domain scores on that row — a single number on the 0–10 scale summarizing mean interference from depression across home chores, school/work, family, and social life. Because rowMeans() uses na.rm = FALSE by default, any respondent missing even one of the four items gets NA on the composite, and the trailing drop_na() then removes them. That leaves a complete-case severity sample of 2,066 adolescents — down from the 2,098 who screened positive, so 1.5% are dropped. This is the sample every Part-2 calculation below uses.

Check the size of what you dropped, every time. Here the fraction removed is small, which limits how far it can move any number — but small does not automatically mean unbiased. What matters as much as the size is whether the missingness is related to the quantities you are studying. Had it come out at 10% or 20%, the honest next question would be who fell out: if the adolescents missing an item differ systematically from those who answered all four, dropping them biases every number that follows. You’d need to say so in the write-up, and the fix would be a principled approach to missing data (i.e., multiple imputation — a topic for PSY 653). The habit to build is checking the fraction and asking the question — not assuming the answer will always be this comfortable.

A note about composite scores

Averaging items assumes the four can reasonably be treated as indicators of one construct, and that equal weighting represents it sensibly. No single statistic settles that: Cronbach’s α speaks to internal consistency, but it does not establish that items measure one thing, nor that equal weighting is right. Construct → operationalization → scoring rule, with a judgment call at every arrow.

On missing items: we required all four items to be present — the complete-case rule. A common, less restrictive alternative averages the available items whenever at least 3 of the 4 were answered (na.rm = TRUE inside a guard that checks the item count) — the “≥ 75% of items” rule you met in M05. Both are defensible; we use the simpler complete-case version here.

Checkpoint 5 · Composite score added

In your notebook — your # Build a composite severity score section now holds this step’s code and output.

Check it worked — Your nsduh_mde_2019 now has a severity_scale column — the row-wise mean of the four domain items on the 0–10 scale — and after drop_na() it holds 2,066 rows.

Step 6 · Describe the distribution empirically

In your notebook — this step fills your # Describe the severity distribution section.

Before fitting any model, look at what you have. What’s the center? How spread out is it? Is it symmetric or skewed?

Get the numeric summary first, then a basic histogram.

Numeric summary

Targets. Pipe nsduh_mde_2019 into a single summarize() that returns, in one row: the sample size, the mean and standard deviation of severity_scale, and the five numbers that describe its shape — minimum, 25th percentile, median, 75th percentile, and maximum. Name the columns whatever reads clearly to you. (quantile() takes the probability as its second argument.)

  • All eight columns go inside one summarize(), separated by commas — the result is a single row.
  • n() gives the sample size and takes no arguments; every other function here takes severity_scale.
  • Center and spread are mean() and sd(); the two ends are min() and max().
  • The middle is median(). For the quartiles, quantile() takes the probability as its second argument: quantile(severity_scale, 0.25) and quantile(severity_scale, 0.75).

The numbers that matter most for normal-model fitting: mean (center), sd (spread), and the min / quartiles / max that describe the shape. Note how close the mean and median sit — that is your first hint about symmetry; Steps 9 and 10 give you two better visual checks.

Your # Setup chunk already loads skimr, and nsduh_mde_2019 |> select(severity_scale) |> skim() returns the same statistics plus a missingness count and a tiny inline histogram. It is the faster move in your notebook, and it’s what the M05 Module used.

We spell the summary out with summarize() on this page because skim() prints two separate tables, and the browser sandbox above renders only one of them — you would see half the output. In RStudio you get both, so use whichever you prefer in your notebook.

Visual summary — a histogram

Targets. Build a histogram of severity_scale from nsduh_mde_2019.

  1. Map severity_scale to x inside aes().
  2. Add geom_histogram(). Set a binwidth deliberately rather than accepting the default — the composite is a mean of four integer items, so it lands on quarter-point values, and your bin width should be a multiple of that. Try one, look at it, and adjust.
  3. Label the axes and give the plot a title that names the variable, the population, and the year.

Map severity_scale to x in aes(). A binwidth of 0.5 groups scores into half-point bins, which matches the quarter-point granularity of a mean of four 0–10 integer items.

Checkpoint 6 · Distribution described

In your notebook — your # Describe the severity distribution section now holds this step’s code and output.

Check it worked — Your notebook shows a numeric summary of severity_scale (mean, SD, quartiles — via summarize() or skim()) and a histogram of the composite score, plus your one-sentence read on its shape.

Write it up — one sentence describing what you see: roughly symmetric or skewed? Is there a single peak, or multiple? Any suspicious clumps at the boundaries?

Step 7 · Fit a normal model — estimate \(\hat{\mu}\) and \(\hat{\sigma}\)

In your notebook — this step fills your # Fit a normal model to severity_scale section.

A normal model has two parameters: the mean \(\mu\) and the standard deviation \(\sigma\). Fitting the model to your data means computing those two numbers from the data and treating them as your best guesses at the underlying \(\mu\) and \(\sigma\). The two estimates are traditionally written \(\hat{\mu}\) (mu-hat) and \(\hat{\sigma}\) (sigma-hat); the hats are the universal statistical notation for “this is an estimate, not the truth.”

\(\hat{\mu}\) is the same number as \(\bar{x}\) from the M06 Module’s parameters-vs-statistics table — the hat notation just emphasizes that we’re using the sample mean as an estimate of the population parameter \(\mu\). Likewise, \(\hat{\sigma} = s\), the sample standard deviation. Same arithmetic, different emphasis: \(\bar{x}\) and \(s\) describe your data; \(\hat{\mu}\) and \(\hat{\sigma}\) use those same numbers as estimators of the underlying normal model’s parameters. Put plainly, we plug in the sample mean and sample SD as our estimates of the model’s \(\mu\) and \(\sigma\).

Targets. Compute the two estimates and store them under the names mu_hat and sd_hat — later steps call them by those names, so match them exactly. Each is one line: pull the severity_scale column out of nsduh_mde_2019 and pipe it into the relevant summary function. Print both when you’re done so you can read \(\hat{\mu}\) and \(\hat{\sigma}\) off the output.

Each line pulls the severity_scale column and pipes into mean() / sd() — e.g. pull(severity_scale) |> mean(). Because you already drop_na()ed in Step 5, there’s no missingness to worry about.

You now have a Normal(\(\hat{\mu}\), \(\hat{\sigma}\)) model of the severity score. Every probability statement you make from here forward uses those two numbers.

Checkpoint 7 · Model parameters estimated

In your notebook — your # Fit a normal model to severity_scale section now holds this step’s code and output.

Check it worked — Your notebook prints two numbers — mu_hat (the sample mean of severity_scale) and sd_hat (its standard deviation) — that together define your fitted Normal(\(\hat{\mu}\), \(\hat{\sigma}\)) model. Every probability in Step 8 uses these two values.

You have a model. Should you trust it yet?

Step 8 is about to use mu_hat and sd_hat to make probability statements — without having looked at whether a normal curve actually describes this variable. In a real analysis that check comes first. If you want to do it in the right order, take the optional fit-checking bonus now, then come back; otherwise carry on, and Step 8 will check the model’s answers against the data instead.

Step 8 · Use the fitted model to make (and check) probability statements

In your notebook — this step fills your # Probability statements from the fitted normal section.

Even if the normal model fits imperfectly, it’s still useful — as long as you know its limits. Now you’ll use pnorm() and qnorm() with your estimated \(\hat{\mu}\) and \(\hat{\sigma}\) to answer two concrete questions, then compare the normal-model prediction to the raw empirical answer from the data itself.

Question 1 — what share of MDE-positive adolescents score above X?

Pick a moderately high severity cutoff on the 0–10 scale, say 7. The normal model says: under our fitted Normal(\(\hat{\mu}\), \(\hat{\sigma}\)), the probability of scoring above 7 is …

Targets. Compute two numbers and compare them.

  1. The model’s answer. Use pnorm() with your mu_hat and sd_hat to get the probability of scoring above 7 under the fitted normal. Store it as p_model. (Watch the direction — pnorm() gives you the lower tail unless you tell it otherwise.)
  2. The data’s answer. Use the emp_cdf function given below to get the observed proportion scoring above 7. Store it as p_empirical.
  3. Print both.

The ecdf() line is given — building it is not what this step is testing. It returns a function: call emp_cdf(7) and you get the proportion of the data at or below 7.

For the model line, q is the cutoff 7, mean/sd are your mu_hat/sd_hat, and lower.tail = FALSE flips pnorm() to “probability above.” For the empirical line, ecdf() gives \(P(X \leq q)\), so 1 - emp_cdf(7) is the matching \(P(X > 7)\) from the raw data.

Compare p_model and p_empirical. Are they close (say, within a couple percentage points of each other)? If so, the normal model is giving you a reasonable answer for this question, even if the fit isn’t perfect overall. If they differ substantially, the normal approximation is poor for this question — and if your goal is to describe this analytic sample, report the empirical proportion rather than the model’s. Note what that does and does not buy you: the empirical proportion is itself a sample statistic, not a population truth. It is the better description of these respondents, which is all either number was ever entitled to claim.

Question 2 — is the model’s 90th percentile actually the 90th percentile?

The normal model claims a specific score is the 90th percentile — meaning 90% of people are at or below that score. You can check that claim directly against the data: under the empirical CDF, what share of the sample actually lives at or below the model’s predicted quantile? If the normal fits, that share should be about 0.90. If it’s not, the gap tells you how far off the model is at that point of the distribution.

Same tool as Question 1 — build the empirical CDF with ecdf(), then query it. In Question 1 we asked “what share is above 7?” → 1 - emp_cdf(7). Here the question is “what share is at or below q_model?” → emp_cdf(q_model) — directly, no subtraction needed.

Targets. Turn the question around — go from a probability to a score, then check the score.

  1. Use qnorm() with your mu_hat and sd_hat to get the score the model calls the 90th percentile. Store it as q_model.
  2. Ask the data what proportion actually falls at or below that score, using the given emp_cdf. Store it as p_empirical_at_q_model.
  3. Print both. If the normal fits well, the second number lands near 0.90.

As in Question 1, the ecdf() line is given.

qnorm()’s p is the probability whose quantile you want — 0.90 — with mean/sd set to mu_hat/sd_hat. Then query the CDF at that score with emp_cdf(q_model).

Two numbers come out. q_model is the score the normal model claims is the 90th percentile. p_empirical_at_q_model is what the data says about that score — the share of observations actually at or below it. If the normal model fits well, the empirical share should be close to 0.90 — but not necessarily exactly 0.90, because we have a finite sample and the composite takes repeated quarter-point values, so small gaps are expected and shouldn’t be over-read.

Checkpoint 8 · Model vs. empirical

In your notebook — your # Probability statements from the fitted normal section now holds this step’s code and output.

Check it worked — Question 1’s two numbers should come out close to each other: the model puts 0.28 of the distribution above 7, and the data put 0.27 there. Question 2’s p_empirical_at_q_model should land near 0.90 — that is the check: the model claims a 90th percentile, and the data say what share actually falls at or below it. If yours is far from 0.90, re-check mu_hat and sd_hat before reading anything into the gap.

Write it up — one or two sentences comparing the two: is the empirical share higher or lower than 0.90, and what does that imply about whether the normal model is placing the 90th percentile in the right spot?

Lab debrief · 5 minutes

Before you do the final render, save your work and look up. We will spend the last five minutes pulling the lab together.

Take one quiet minute to skim the questions below, then we’ll discuss as a group.

Lab debrief · what did we learn by doing?

  1. The sticking point. What was the single hardest moment in today’s lab — the place where you got stuck the longest? Was it reading a probability off the cross-tab, or judging whether the normal model actually fit? What helped you move forward, or what would have helped?

  2. Reading a 2×2 table as six probabilities. Every marginal, joint, union, and conditional probability lived somewhere in the rows, columns, or interior cells of your cross-tab. Which of the six was hardest to see in the table, and which cell or margin did you have to divide by to get it?

  3. The direction of the bar matters. \(P(\text{SUD} \mid \text{MDE})\) and \(P(\text{MDE} \mid \text{SUD})\) came out to different numbers even though they share the same joint in the numerator. Which one would a clinician screening for co-occurring conditions care about more, and why? Where else in your own research would getting the direction backwards change the claim?

  4. Looking before believing. Part 2 had you fit a normal model and then compare what it predicts against what the data actually show. How different were the model-based and empirical answers in Step 8 — and what would you report if you had to write this up? (The optional fit-checking bonus adds the two visual checks that belong between those steps.)

Final render and submit · at home

  1. Add your name to the author: field at the top.
  2. Do a final render. Click Render or press Cmd/Ctrl + Shift + K.
  3. Open the rendered file — it renders next to the .qmd in programs/.
  4. Read it end to end — open the HTML in a browser and read it as if you were a stranger.
  5. Submit it to Canvas under “Lab 6 — Probability Practice with NSDUH.”
  6. Commit and push. In GitHub Desktop, commit m06_lab.qmd with a one-line summary of what you added, then Push origin.

Double check

Before you leave today:

Every piece of prose the lab asked you to write — three sections in all. A notebook of outputs with no writing is a pile of results, not an analysis:

What you just did, in research terms

You took a real NSDUH data file and used probability in two different ways, with the same vocabulary behind both. In Part 1 you read the co-occurrence of depression and substance use straight off a 2×2 table, the exact bookkeeping behind any sentence that begins “adolescents with MDE were more likely to…” In Part 2 you fit a normal model to a severity score, checked it before trusting it, and used it to make tail-probability statements you then sanity-checked against the raw data. That “estimate → check → use with appropriate caution” arc is the whole engine of applied statistics.

What you have not done yet is statistical inference. Both halves of today described this sample: Part 1’s probabilities are unweighted proportions in the complete-case analytic subset, and Part 2 fitted a model to an observed distribution and compared it against that same data. Nowhere did you quantify how much any of those numbers would move if the sampling process were repeated — and that question, not the arithmetic, is what inference means.

M07 asks it: how precisely does one sample’s \(\bar{x}\) pin down the population value it estimates? Once a statistic is treated as a random variable with its own sampling distribution, every probability tool from this Module can be pointed at that distribution — and probability becomes the machinery of inference. You’ll build your first confidence intervals there, by bootstrap and then with the t-distribution.

More practice (optional)

After class, or before an exam. Nothing here goes into your lab notebook, and none of it is graded.

Three standalone pages you can run any time — each drills one skill from today on different data. Two extend Part 1, one extends Part 2.

Every optional activity in the course is also listed on one page: Optional Activities.