Reproducing the Illusion of Predictability

Lab · Module 9 · Wed Oct 14

Welcome to the Module 9 lab

Hypothesis-testing week ends with real published data. Today you run a computational reproduction: take a paper’s own data, re-run its analysis, and see whether the reported numbers come back out. (The other credibility check is a direct replication — collect new data under the same design and see whether the finding recurs. Reproduction asks “can this result be recovered from these data?”; replication asks “does it recur in new data?” — and they fail for different reasons.) Reproduction is the one you can do in a class session, and it is also the shape of Project 2: there you find the paper, procure the data, and defend the analytic decisions. Today all three are handed to you.

The paper is An illusion of predictability in scientific results: Even experts confuse inferential uncertainty and outcome variability (Zhang et al., 2023, PNAS) — open access, with a copy in the course readings folder as zhang_2023.pdf and the authors’ data and code at github.com/jhofman/illusion-of-predictability, where our three course datasets come from. It is co-authored by the same Dr. Jake Hofman behind the boulder-sliding study, and you met its headline numbers in the M07 Module. Same question — do experts confuse a chart’s inferential uncertainty with its outcome variability? — now asked of clinicians, data scientists, and faculty.

Note

Two papers this week — keep them straight. The Module and pre-study use Hofman et al. (2020), the boulder study, with Experiment 1 and Experiment 2. Today’s is Zhang et al. (2023), which reports Study 1, Study 2, and Study 3 on three expert populations — different paper, different participants, same underlying question.

What you’ll leave with

  • A reproduced Welch’s t-test and Cohen’s d for one expert-audience study, run with infer and effectsize
  • A picture of your study’s results showing error bars and individual points — built to the standard the paper argues for
  • An estimation plot putting the effect and its uncertainty on their own axis — a figure you can reuse for Project 2
  • A completed APA result line, your numbers next to the paper’s, and a reproduction verdict
  • A shared 3-row reproduction checklist for the paper’s focal result across all three studies, built with the class

This is a jigsaw lab: your group reproduces one study, then we reconvene so all three fit together. Each analysis step has ✍️ Your Code, 💡 Hint, and 👀 Spoiler tabs; as in M06–M08, the ✍️ box starts empty and a Targets list above it says what to build. One figure step hands you the plot’s scales and labels and asks only for the layers that carry the statistics; the estimation plot is given to run and read. Stuck for more than a few minutes? Open 💡, then 👀, and keep moving — your group is waiting on your numbers for the share-back.

Before you start — set up your notebook

As in M07 and M08, there is no per-lab notebook template: copy lab_template.qmd and build the structure yourself — Step 0 scaffolds it in a few minutes. This page’s sandboxes are the group scratchpad; your own notebook is the individually-submitted deliverable. The three illusion datasets are already in your project’s data/ folder, each with a codebook in documentation/. Steps 1–4 and the debrief happen in class; the write-up — the interpretation prose, the final render, the Canvas submission — is finished at home, and each step ends with an In your notebook note saying exactly what to add.

Carrying code into your notebook — the same rules as M06

The code boxes on this page run in the browser; your notebook runs on your copy of R. So anything you carry across has to sit inside an R chunk you insert yourself, with a label on a #| line. If code lands outside a chunk it renders as plain text and never runs. And your notebook runs top to bottom in a fresh session every time you render, so paste the chunks in the page’s order.

The M06 lab’s Step 0 walks through the keystrokes for inserting and labelling a chunk, and what goes wrong when code lands outside one — worth a look if any of that is hazy.

Step 0 · Get set up

Fourth time through — a few minutes. (The M06 lab’s Step 0 has the line-by-line explanations.)

  1. GitHub Desktop first. Select PSY652_project, click Fetch origin, then Pull origin if it appears — the pull half of the loop you’ll close at the end of the lab.

  2. Open RStudio via PSY652_project.Rproj, and open programs/lab_template.qmd.

  3. File → Save As… a copy in programs/ named m09_lab.qmd. The template opens with three stock sections — # Setup, # Import Data (which reads nhanes purely as a worked example), and # Glimpse. Delete all three: you replace them with the setup and import chunks below, and two chunks with the same label stop a render.

  4. Set the YAML — change the title, and bump toc-depth to 5 so the report’s ## subsections reach the table of contents:

    title: "Reproducing Zhang et al. (2023): The Illusion of Predictability"
toc-depth: 5

  5. Give the notebook this section skeleton — you fill it as the lab proceeds:

    # Introduction
# Data
# Results
## Picture of the data
## The effect, drawn
## Welch's t-test and effect size
## APA result and reproduction verdict
# Discussion

  6. Add the setup chunk — packages only, right after the YAML, with a label so you can find it:

```{r}
#| label: setup

library(tidyverse)
library(infer)
library(effectsize)
library(gt)
```
  1. Import the data in its own chunk under # Data. Packages and data stay in separate chunks on purpose: an error in setup means a missing package; an error here means a wrong path.
```{r}
#| label: import-data

# Medical providers -- we reproduce the blood-pressure scenario
medical_bp <- read_rds(here::here("data", "illusion_medical.Rds")) |>
  filter(scenario == "Blood pressure scenario") |>
  mutate(condition = fct_relevel(condition, "Saw SDs first", "Saw SEs first"))

# Data scientists
datasci_analysis <- read_rds(here::here("data", "illusion_datasci.Rds")) |>
  mutate(condition = fct_relevel(condition, "SE + Points", "SE Only"))

# Faculty
faculty_analysis <- read_rds(here::here("data", "illusion_faculty.Rds")) |>
  mutate(condition = fct_relevel(condition, "SE + Points", "SE Only"))
```

All three files load even though your report analyzes only your group’s — that keeps the share-back comparisons reproducible. You may drop the two you didn’t use once the share-back is done.

Step 1 · The paper in three minutes

Zhang et al. (2023) · An illusion of predictability in scientific results

Citation. Zhang, S., Heck, P. R., Meyer, M. N., Chabris, C. F., Goldstein, D. G., & Hofman, J. M. (2023). An illusion of predictability in scientific results: Even experts confuse inferential uncertainty and outcome variability. Proceedings of the National Academy of Sciences, 120(33), e2302491120.

The central claim. When a display emphasizes inferential uncertainty — here, standard-error bars — while making outcome variability hard to see, viewers can substantially overestimate effect sizes and how predictable individual outcomes are, even when the viewers are experts. In Study 1, showing SDs reduced the large overestimates SE displays produced, though those estimates tended to undershoot the true effects. In Studies 2 and 3, showing the individual outcomes alongside the SEs produced well-calibrated judgments. The authors therefore recommend showing outcome variability alongside statistical estimates when possible.

One focal question, all three studies. Participants saw the results of a hypothetical medication trial (Study 1) or an extended abstract based on a published study of violent video games (Studies 2 and 3), then estimated the probability of superiority (PSup) on a 50–100 scale: pick one person from the treatment group and one from the comparison group — how likely is the treated person to have the better outcome? 50 means no advantage either way; 100 means a sure thing. (PSup was not the only outcome collected, but it is the one measured the same way in all three studies, and it is the one you reproduce.)

One visual feature changed.

  • Study 1 — Medical Providers. Group means with SE bars versus the same means with SD bars. One twist: every clinician eventually saw both figures — what was randomized was which came first, and the focal comparison uses each clinician’s first estimate. That is why the conditions are named Saw SEs first and Saw SDs first, and why comparing them is still a clean randomized two-group comparison. (Group A analyzes the blood-pressure scenario, the Study 1 comparison we reproduce; the paper reports the COVID-19 scenario alongside it.)
  • Studies 2 and 3 — Data Scientists and Faculty. Straightforward between-subjects designs: SE Only versus SE + Points — the same means and error bars plus the individual observations.

Two independently randomized groups, one continuous outcome, a question about two population means — the M09 decision map points to a two-sample Welch’s t-test, applied three times. If the authors are right, the groups that could not see outcome variability should report higher PSup.

The three studies you will reproduce:

Study Population Visualization contrast Course-file N Published focal test
Study 1—Medical Providers clinicians SE bars versus SD bars, blood-pressure medication scenario 75 t(65.1) = 6.83, p < .001
Study 2—Data Scientists professional data scientists at a software company SE-only bars versus SE bars plus individual points 175 t(159.4) = 6.34, p < .001
Study 3—Faculty tenure-track academics SE-only bars versus SE bars plus individual points 368 t(363.9) = 4.52, p < .001

The N column is computed from the course files you are about to load. The final column is quoted from the published paper. Those are different kinds of numbers: one comes from the dataset in front of you; the other is the published result you are attempting to reproduce.

How close should your numbers come?

Our course datasets come from the authors’ public repository, analyzed the way the paper reports: R’s default two-sided Welch t-test. Participants with missing PSup estimates were removed when the course files were built, so the N in each file is the N analyzed today.

  • Study 1: expect t(65.1) = 6.83 — matching exactly.
  • Study 3: expect t(363.9) = 4.52 — matching exactly.
  • Study 2: expect about t(159.2) = 6.41, against the paper’s t(159.4) = 6.34. The published analysis applied a preregistered background-survey exclusion our course file doesn’t implement — and in the authors’ own code that exclusion removes exactly one participant, who is the entire gap. A near-match with a documented reason is a successful reproduction; Group B states the difference in its verdict. (“Why one, when the Methods say two?” — the box below walks through it, and it is the clearest example in this course of why code is part of the scientific record.)

Published statistics depend not only on the final test but on the decisions that built the analytic sample. “Our statistic was 6.41 rather than 6.34 because our file retains one participant the published analysis excluded” is more scientifically informative than a match you cannot explain.

When the prose and the code disagree

Where this comes from: the article itself establishes one side of the discrepancy — that two participants were removed. Everything else in this box comes from inspecting and running the authors’ released analysis code.

Here is what a computational reproduction check is actually for, and Study 2 gives us a remarkably clean example. In short: the paper’s Methods say two participants were removed, but running the authors’ own published script removes one — and the reason is a silent R type coercion that no amount of re-reading the Methods could reveal. That single participant is the whole 6.41-vs-6.34 gap.

The full walk-through is below. It is optional — Group B needs only the one-participant fact above — but it is the clearest example in this course of why code is part of the scientific record.

The paper’s Methods report that the authors removed two participants. But the published analysis script builds the exclusion object like this — quoted from analysis/shared_analysis.R, de-indented (it sits inside a function) and still in the older %>% pipe rather than the native |> this course uses:

assignments_to_drop1 = background_ds %>%
  filter((has_any_stats_training == 0) | (rctComfort == 1)) %>%
  select(assignmentId)

assignments_to_drop = c(
  background_with_number_done %>%
    filter(number_done == 0) %>% select(assignmentId),
  assignments_to_drop1)

(background_with_number_done is the background survey file with one column added: a per-person tally of how many experience categories they reported, so number_done == 0 finds anyone who reported none.)

Between them the two filters identify three unique participant IDs. The first identifies one participant; assignments_to_drop1 identifies three, including that same participant.

The problem is the use of c(). Each input is a one-column data frame, and c() does not stack their rows into one vector of IDs. Instead it creates a list with two elements:

  1. one element containing a single ID;
  2. one element containing a vector of three IDs.

That structure produces two subtle consequences.

First, the script counts the number excluded with:

background_numbers["num_dropped"] = length(assignments_to_drop)

The length of this object is 2 — because it holds two list elements, not because it holds two participant IDs. That count is saved to results/das_results.Rdata, the file the manuscript reads its numbers from, which appears to explain the figure reported in the paper’s Methods.

Second, the analysis later filters participants with:

filter(!(assignmentId %in% assignments_to_drop))

The first list element holds one bare ID, so that participant is matched and removed. The second holds three IDs bundled as a single vector, which coerces to the literal text c("id1", "id2", "id3") and matches no individual assignmentId. The published code therefore removes one participant, not the three the filters identify.

The numbers confirm it:

Analytic sample n Welch test
No participants removed — our course file 175 t(159.2) = 6.41
Published code — one participant removed 174 t(159.4) = 6.34
All three flagged participants removed 172 t(155.3) = 6.26

The middle row is the statistic printed in the paper. The entire difference between our result and the published one is therefore a single participant — and a quiet data-structure error that becomes visible only when the code is executed and its intermediate objects are inspected.

The intended construction would be something like:

assignments_to_drop <- bind_rows(
  background_with_number_done |>
    filter(number_done == 0) |>
    select(assignmentId),
  assignments_to_drop1
) |>
  distinct(assignmentId) |>
  pull(assignmentId)

Two things to take from this:

  1. The discrepancy does not threaten the study’s conclusion. Whether zero, one, or all three flagged participants are removed, the estimated difference stays large and the evidence against equal group means stays overwhelming. The preprocessing error moves the reported statistic slightly; it does not move the finding.
  2. Code is part of the scientific record. A Methods section describes the analysis the researchers intended to run; executable code reveals the analysis the computer actually ran. A strong reproduction does not merely ask whether the final number matches — it traces any difference back through the data-processing decisions that produced it.

When you complete Project 2, write your code as though another researcher will run it — because reproducible science assumes that eventually someone will.

Your assignment

We’re splitting into three groups, one per study:

  • Group A → Medical Providers (Study 1)
  • Group B → Data Scientists (Study 2)
  • Group C → Faculty (Study 3)

In Step 3’s section a you’ll pick your study’s tab and name its two conditions; sections b–f are then the same pipeline for everyone, run on your own study’s data. We reconvene in Step 4 to share back.

Checkpoint 1 · You can frame the reproduction

Before you touch code, you can say in one sentence each: (a) the paper’s central claim — even experts confuse a chart’s inferential uncertainty with its outcome variability; (b) what your assigned study manipulated; and (c) that every study’s focal test is a two-sample Welch’s t on a 50–100 probability-of-superiority estimate. You know which group — A, B, or C — is yours.

In your notebook — your Introduction borrows straight from this step: two or three sentences giving the paper’s question and what your group’s study manipulated.

Step 2 · Whole-class setup

We load all three datasets up front so everyone shares one starting state. Each .Rds is the tidied analytic dataset for one study, prepared from the paper’s public-repo CSVs.

One decision before anything else: which condition is the baseline. A factor’s level order silently decides which group R treats as the reference — and therefore the sign of every difference you report. In all three studies we put the group that saw outcome variability first, so a positive difference always means the error-bars-only group estimated higher: the direction of the sentence you’ll write.

Dataset Baseline (saw outcome variability) Order as saved
illusion_medical Saw SDs first Saw SEs first, Saw SDs first needs flipping
illusion_datasci SE + Points SE + Points, SE Only already first
illusion_faculty SE + Points SE + Points, SE Only already first

Only the medical file actually needs reordering, but the chunk calls fct_relevel() on all three anyway — depending on whatever order a file happened to be saved in is fragile, and stating the order you want is how you make sign-flips impossible. The same chunk filters the medical file to the blood-pressure scenario, the Study 1 comparison we reproduce.

Two surfaces, two file paths. In the sandbox below, the data are mounted at shared_data/, so the path is a plain string and there is no here package. In your own notebook the same files live in your project’s data/ folder, so you write here::here("data", "illusion_medical.Rds") — exactly as in the Step 0 import chunk you already pasted. Same files, different location, so the path differs. If you paste one version into the other, R will tell you it can’t find the file.

illusion_medical / illusion_datasci / illusion_faculty · 163 / 175 / 368 observations · shared schema · Zhang et al. (2023)

Three tidied analytic datasets — one per expert population — prepared from the paper’s public-repo CSVs (the medical file spans both scenarios; our focal reproduction uses the 75 blood-pressure rows). They share a single analytic schema:

  • participant_id character — anonymized worker ID
  • condition factor — the chart the participant was shown (levels differ per study; Step 3’s section a names yours)
  • psup_estimate numeric — the participant’s estimate of probability-of-superiority, on a 50–100 scale
  • scenario factor — medical only: which medication scenario the participant was assigned (blood-pressure vs COVID-19)
  • study character — constant within each dataset (medical / datasci / faculty)

Before you test: what the t-test needs

psup_estimate is bounded — it cannot leave 50–100 — and piles against those edges, so the individual scores are not Normal. That’s fine: what the test needs is for the sampling distribution of the mean difference to behave, and at these per-group sizes — from 34 in Study 1’s smaller arm to 199 in Study 3’s larger one — M07’s Central Limit Theorem makes the procedure reasonably robust, especially because the bounded scale rules out arbitrarily extreme observations that could drag a mean around. Independence comes from the design: one estimate per person, each randomly assigned to one condition. Welch’s version also declines to assume the two groups share a variance, which is why its df comes out fractional, like 65.1. The habit worth building: not “is my outcome Normal?” but “is the mean well estimated at this sample size, and are the observations independent?”

Checkpoint 2 · Data loaded and understood

All three datasets load, and your glimpse() shows the shared schema — condition and psup_estimate are the two columns every test needs. You can say what one row represents — for Study 1 that is a clinician’s first estimate; for Studies 2 and 3 it is each expert’s focal Part 1 estimate.

In your notebook — this load already lives in the import chunk you added under your Data heading in Step 0 (note the path differs there: here::here("data", ...), not the sandbox’s shared_data/). Add a sentence to that Data section naming your group’s dataset and what one row represents; the glimpse() is a check you don’t need to keep.

Step 3 · Reproduce your study

How to use this part

Pick your study in a — that tab sets my_data, your two condition names, and your study’s title. Sections b → f are then the same pipeline for everyone, run on your my_data, so every Spoiler prints your study’s numbers. Work in order; if you’ve been stuck more than a few minutes, open 💡, then 👀, and keep moving.

a · Pick your study

Everything on this page runs in one shared R session, so opening another group’s tab and running its code replaces my_data for the rest of the page. To put it back, re-run your own tab and every chunk after it, section b’s summaries included; the estimation plot stops with a message if they don’t match. (The Analyzing: line each tab prints, and your study’s name on every plot title, are the tell.)

This section holds the one genuinely analytic decision of the pipeline: naming which condition is the baseline (ref_group) and which is compared against it (comp_group). That choice fixes the sign of the difference, the t, and the d — use the exact labels count() prints, not memory.

Your dataset is medical_bp. Step 2 filtered it to the blood-pressure scenario (the Study 1 comparison we reproduce) and put “Saw SDs first” — the group that saw outcome variability — first, making it the baseline. Each clinician eventually saw both chart types; the paper compares their first estimates, which is why the conditions are named “Saw SEs first” and “Saw SDs first”.

The baseline is the group that saw outcome variability — SD bars or individual points, depending on your study. Both strings must match the printed factor levels exactly, including capitalization and spacing.

Your dataset is datasci_analysis. The conditions are “SE Only” (just error bars) and “SE + Points” (error bars plus individual data points). The paper’s claim: the SE-only group will overestimate the probability of superiority because they can’t see the spread of individual outcomes. A heads-up before you run the test. Your tobs will land near 6.41 — close to, but not exactly, the paper’s reported t of 6.34. That’s expected: the authors applied a preregistered exclusion keyed to a background file our course dataset doesn’t ship, and the entire gap comes down to one participant. A near-match with a documented reason is still a successful reproduction check — and this is exactly the kind of divergence your Project 2 write-up should notice and explain rather than hide. (The “When the prose and the code disagree” box back in Step 1 has the full story, and it is worth reading before you write your verdict.)

The baseline is the group that saw outcome variability — SD bars or individual points, depending on your study. Both strings must match the printed factor levels exactly, including capitalization and spacing.

Your dataset is faculty_analysis. Same conditions as Study 2 — “SE Only” vs “SE + Points” — but the sample is tenure-track academics drawn from psychology, sociology, physics, biology, business, and computer science.

A documentation wrinkle in this study too, worth one sentence in your verdict. Your test will reproduce the paper’s t(363.9) = 4.52 exactly, using all 368 responses. The article also reports that 63 participants were excluded. Those two statements cannot both describe the same analysis: Welch’s degrees of freedom can never exceed N − 2, so a 305-person sample could not produce df = 363.9. We reproduce the published statistic as reported; how the stated exclusion count relates to that particular test isn’t clear from the article. Note it and move on — it does not change the finding.

The baseline is the group that saw outcome variability — SD bars or individual points, depending on your study. Both strings must match the printed factor levels exactly, including capitalization and spacing.

From here on every step runs on my_data — your numbers will differ from the other groups’, and that is the point of the share-back.

In your notebook — the first chunk under # Results is your group’s section-a chunk, copied whole. It creates the four objects every later chunk uses — my_data, study_title, ref_group, and comp_group — so copying only the my_data line leaves the next chunks with nothing to work on. Copy just your own group’s tab: your notebook analyzes one study.

b · Descriptives, then the picture

First, the numbers your APA line will need.

Targets. Build one row per condition, each carrying a mean and its interval.

  1. Start from my_data, group_by() condition, and summarize().
  2. Six columns: the cell size, the mean of psup_estimate, its SD, its standard error, and the two 95% bounds.
  3. The SE of a mean is the SD over the square root of n — M07’s quantity, now in service of a comparison. For the bounds, the multiplier is t* from qt(); the Hint names the call.
  4. Name the columns n, M, SD, SE, lo, hi — the figure refers to them by those names. Set .groups = "drop", assign to group_descriptives, and print it.

The standard error of a mean is the SD divided by the square root of the group size — you have both as columns already, so SD / sqrt(n).

For the bounds, the multiplier is not 1.96. That value is the large-sample shortcut; the model-based multiplier comes from the t-distribution with this group’s degrees of freedom, which qt() gives you: qt(0.975, n - 1). Note 0.975, not 0.95 — a two-sided 95% interval leaves 2.5% in each tail.

Column What it is What it’s for
n how many participants in that condition — a group’s n that group’s SE and its interval’s df. Add the two groups’ n for the total N in your APA line
M the group mean the M in your APA line
SD how much individual participants differ from each other the SD in your APA line
SE how precisely that group’s mean is estimated — SD / sqrt(n) the width of that group’s interval. Welch’s test later combines both groups’ uncertainty into the standard error of the difference — neither group’s SE is itself the test’s denominator
lo, hi the 95% confidence interval on that group’s mean the error bars in the chart below

The distinction that matters: SD is about people, SE is about the estimate. They answer different questions, and only one of them shrinks as you collect more data. If your group hesitated anywhere, it is almost always here.

Now plot your study’s results. This figure shows what your participants estimated — it is not a recreation of the stimulus figures Zhang et al. showed them. And its design follows the paper’s own recommendation: a results figure with error bars alone would commit, in your write-up, precisely the error the paper documents. So it needs three layers — every individual observation (jittered), each group’s mean, and each group’s 95% interval from the lo/hi you just computed.

One construction detail: layers 2 and 3 come from group_descriptives (two rows) while layer 1 comes from my_data (one row per person). A single ggplot() can draw from two data frames — pass data = inside the geom, with inherit.aes = FALSE so that layer ignores the plot-level aes().

Targets. The scales, zoom, and labels are given — add the two layers that carry the statistics: the jittered raw observations, and the error bars mapped to your bound columns.

For the individual points you want the jittered scatter from M03 — geom_jitter(), which nudges points sideways so ties don’t hide each other. (Plain geom_point() would stack them into a single vertical line.)

The error-bar layer reads from the summary table you built above, and geom_errorbar() wants a bottom and a top: ymin and ymax, which are exactly the two bound columns in that table.

Read it before moving on: the cloud is the people, the dot is the group mean, and the bar is how well you know that mean.

A note on what we’re plotting, versus what the paper plots. Zhang et al.’s figures show the mean with one standard error above and below it. We draw 95% confidence intervals instead, because they answer the question a reader actually has — what range of population means is this sample consistent with? — and because it keeps this chart consistent with the interval you build in section c. So your figure is deliberately not a pixel match for theirs. Say so if you reproduce it in your write-up.

c · Draw the effect

The plot above shows the two groups, but a reader is still left to infer the difference between them — and how precisely it is known — from two separate summaries. An estimation plot (Ho et al., 2019) makes that contrast explicit while keeping the individual variability in view, a design that fits the paper’s broader lesson that inferential summaries and outcome variability shouldn’t be visually conflated. It puts three things in one frame:

  1. Every observation, so nobody takes a summary on trust.
  2. Each group’s mean with its 95% interval, in the outcome’s own units.
  3. The difference itself, on its own axis at the right, with its own uncertainty — equal visual weight, instead of something the reader computes by eye.

The right-hand axis is the left-hand ruler relabeled so zero sits at the baseline mean (sec_axis(~ . - baseline)), and the figure deliberately mixes two kinds of 95% interval — t-intervals on the group means, a percentile bootstrap on the difference. The Going further box unpacks both.

This estimation plot is a course addition, not part of the paper’s focal reproduction: Zhang et al.’s published focal result is the Welch test in section d. (Where the authors bootstrapped elsewhere in the paper, they used 10,000 resamples and a different interval method, the reverse percentile.) The same holds for Cohen’s d and the CLES in section e — useful additions, not reproduction targets.

“Wait — didn’t M03 say two y-axes are bad?” It did, and the rule is a good one. M03’s trap list classes dual y-axes as bad, on the grounds that putting two variables on different y-axes lets you manufacture any correlation you like by adjusting the scales. (The critique is an old and well-argued one — see Few, 2008 — and it is why ggplot2 makes a genuine second axis so awkward to build.)

Read that reason carefully and you’ll see this figure isn’t the thing being warned about. The problem case has two different variables and two independently chosen scales — and it’s that second freedom that does the damage, because sliding one scale against the other can make an association appear, vanish, or reverse.

Here there is one variable (probability of superiority) on one scale. The right-hand axis is not a second measurement; it is the same ruler, relabeled so that zero falls at the baseline group’s mean — which is exactly what sec_axis(~ . - baseline) says: take this axis, subtract a constant. Nothing is free to be tuned, so nothing can be manufactured. Slide the data and both axes move together.

The transferable lesson is about how to hold a design rule. “Never use two y-axes” is a shortcut for “don’t give yourself a free parameter that can invent a relationship.” When you meet a figure that breaks the surface form of a rule, check whether it also breaks the reason. Sometimes it doesn’t — and knowing why is what separates following conventions from understanding them.

Two kinds of interval in one figure. This figure deliberately mixes two procedures, and you should be able to say which is which:

  • The group means on the left carry parametric t-intervals — the qt(0.975, n - 1) bounds you computed in section b.
  • The difference on the right carries a percentile bootstrap interval — resampled, nonparametric, straight out of M07. (Nonparametric means it assumes no distributional shape — but it does still assume your observations are suitable units to resample, which is the independence you get from the design.)

Both are 95% intervals, and here they will tell the same story. They are not the same procedure, though, and they are not guaranteed to agree exactly — the bootstrap interval on the difference need not match the Welch interval your t-test reports in section d. Naming which interval came from which method is part of describing a figure honestly.

The bootstrap, aimed at a difference

The same bootstrap you learned in M07 works here. Back then, each resample produced one mean; here each resample produces a difference between two means. The resampling steps are identical — only the statistic you calculate in summarize() changes. The middle 95% of the 2,000 differences is the percentile bootstrap interval for the difference.

We resample whole participant rows, so each person’s PSup estimate travels with their condition. Because we resample the complete dataset rather than each condition separately, the number of participants in each condition wobbles a little from resample to resample; a bootstrap that resampled within each condition would hold those group sizes fixed. Both give nearly the same interval here, and we use this version because it is exactly the bootstrap workflow you already know.

Notice what the bootstrap does not do: it never imposes a null. The resampled differences are centered on what your data actually show, which is what you want when the question is how big the difference is and how precisely it is pinned down. That is the difference between estimating (this figure) and testing (the t-test in section d), and it is worth remembering in general.

The figure itself is given — run it, then read it:

Starting point

You have your two groups’ summaries from section b and your two condition names from section a. You want one figure that shows every participant, each group’s mean with its interval, and the difference between the groups with its own uncertainty.

What the code does, in four steps

  1. Positions and the gap. as.numeric(condition) turns your two conditions into x positions 1 (baseline) and 2 (comparison), leaving position 3 free for the difference. baseline is the baseline group’s mean, and diff_obs is comparison minus baseline.
  2. The bootstrap — M07’s, with a new statistic. The resampling steps are the ones from the M07 pre-study: rep_slice_sample() draws 2,000 resamples of your whole dataset with replacement (prop = 1 makes each resample the same size as your data; a bootstrap needs that, so use either prop = 1 or n = the sample size, because smaller resamples would make the interval too wide), and group_by(replicate) keeps each resample separate. Only the statistic in summarize() changes: instead of one mean, each resample’s comparison mean minus its baseline mean. Read psup_estimate[condition == comp_group] as “the PSup estimates of the people in the comparison condition.” quantile() then keeps the middle 95% of the 2,000 differences.
  3. The sideways curve. density() smooths the 2,000 differences into a curve, the same smoothing a density plot does. The next lines lift that curve so its zero sits at the baseline mean, and lay it on its side at x = 3.
  4. The plot, bottom layer to top. The jittered participants; a dashed line at the baseline mean; the sideways bootstrap curve; each group’s interval and mean; and at x = 3, the difference’s interval and the difference itself. The last scale line adds the right-hand axis.

Key functions

Function What it does
rep_slice_sample(prop = 1, replace = TRUE, reps = 2000) M07’s bootstrap: 2,000 resamples, each the size of your data, drawn with replacement.
group_by(replicate) |> summarize(stat = ...) One statistic per resample. In M07 it was a mean; here it is a difference of two means.
x[condition == comp_group] Keeps only the values of x for people in the comparison condition.
quantile(stat, 0.025), quantile(stat, 0.975) The cutoffs of the middle 95%: the percentile bootstrap interval.
density() Smooths a pile of numbers into a curve.
geom_ribbon(orientation = “y”) Fills an area sideways: at each height, from x = 3 out to the curve.
annotate(“linerange”, …) Draws one vertical bar at a fixed position: the difference’s interval.
sec_axis(~ . - baseline) A second y-axis: the same scale, shifted so zero sits at the baseline mean.

How to read output

A single figure, plus two printed lines with the difference and its interval. The box just below walks through the figure piece by piece.

Print the difference and its interval — the raw effect and its uncertainty, which lead your APA line:

Reading an estimation plot — a figure worth keeping for Project 2

Read it left to right:

  • The dots are every participant’s estimate, by condition, on the left axis (PSup, 50–100). They show how much individuals vary.
  • The black point and bar on each group are its mean and 95% t-interval from section b — how precisely each group’s mean is known.
  • The dashed line marks the baseline group’s mean. The right-hand axis is the same ruler relabelled so that line reads zero, so anything on the right reads as a difference from the baseline.
  • At the far right, the blue dot is the estimated difference (comparison minus baseline), the vertical bar is its 95% bootstrap interval, and the shaded curve is the 2,000 bootstrap differences, turned on their side.

One picture answers three questions in the outcome’s own units: how much do people vary, how well is each mean known, and how big is the difference — and how precisely is that known? If the difference’s interval clearly excludes zero, you have strong reason to expect the t-test in section d to reject \(H_0\). They are different procedures, so near zero they can part company; here they tell the same story.

For Project 2: any two-group comparison you reproduce can be shown this way. The chunk runs on my_data, ref_group, comp_group, and study_title, so pointing it at your own data means changing those four objects and the outcome’s column name.

Talk it through: read the figure before you test it

Nobody runs the next chunk until the group can answer, from the picture alone:

  • Point at the effect. Which mark on this figure is the estimated difference? Which one is the uncertainty in it?
  • Someone claims the two groups “barely differ because the dots overlap so much.” Point to what they are confusing.

The effect is the colored dot on the right-hand axis; the vertical bar through it is its 95% interval, and the shaded curve behind shows which values the resampling supported most often. The group means on the left are inputs to that difference, not the effect itself.

The confusion is between individual spread and precision of a mean. Overlapping dots say individuals vary a lot. The interval on the difference says the average gap is pinned down well. Both are true at once — the same point the CLES makes in the Going further box under section e.

d · The Welch’s t-test

Hypotheses first. In symbols, with your two groups:

\[H_0: \mu_{\text{comp}} - \mu_{\text{ref}} = 0 \qquad\qquad H_a: \mu_{\text{comp}} - \mu_{\text{ref}} \neq 0\]

Two cautions. The parameters are population means for the two conditions — not any individual’s estimate. And the paper predicts a direction (the group that could not see outcome variability estimates higher) but tests two-sided: predict a direction, test two-sided, and let the estimate carry the sign. In your notebook, write both hypotheses with your study’s actual condition labels.

Targets. Run the focal comparison.

  1. Pipe my_data into t_test().
  2. formula: outcome ~ grouping variable — the outcome is psup_estimate.
  3. order uses the two names from section a, comparison first, so the difference is comparison minus baseline — the same direction as your estimation plot.
  4. Two-sided alternative; assign to test_result and print it.

t_test() from infer runs Welch’s version by default — no equal-variance assumption; the M09 Module explains why that is the right default.

The formula puts the outcome on the left and the grouping variable on the right — the two column names from glimpse(). For order, you stored the right strings as comp_group and ref_group in section a; pass them in that order. The test is non-directional, so alternative is "two-sided".

Read the output in the order the Module taught. estimate is the raw mean difference in percentage points — the number your APA line leads with. lower_ci/upper_ci bound it. statistic is tobs, t_df its degrees of freedom, p_value the two-sided tail area.

Notice that t_df is not a whole number. That is Welch’s degrees of freedom: because the two groups are allowed to have different variances, the df is a weighted blend of each group’s contribution to the standard error rather than a simple count. Compare statistic and t_df with your study’s row in Step 1’s table.

Talk it through: one question

  • Suppose your study had only N = 20 participants in total. Which of estimate, statistic, and p_value would change most? Which would change least?

With N = 20: shrinking the sample does not change the population quantity you are estimating — but be careful, that is not the same as saying your observed estimate would stay put. A 20-person estimate is far less stable; it could land well away from the current value depending on which 20 people you happened to draw. What changes systematically is the standard error: it grows, so |statistic| generally shrinks and p_value generally rises. That asymmetry — the target holding still while the evidence about it weakens — is the difference between an effect and the evidence for it.

e · Effect size (Cohen’s d)

The t-test says the difference is hard to attribute to sampling variation. Cohen’s d says how big it is, in standard-deviation units.

Why pooled_sd = FALSE: Cohen’s d needs a standardizer, and the default pooled SD assumes the two populations share one variance — exactly the assumption Welch’s test declines to make. pooled_sd = FALSE standardizes by \(\sqrt{(s_1^2 + s_2^2)/2}\) — the square root of the average of the two group variances, the Module’s formula — making no equal-variance claim, and it is the package’s own recommendation to accompany a Welch test. The two numbers answer different questions: Welch’s t measures the difference against its sampling uncertainty; d measures it against how much individuals differ. One reporting habit follows: name your standardizer — “d using the unpooled SD” is a complete description; a bare d is not.

A good first addition when you assemble your notebook at home. Cohen’s d is in standard-deviation units — precise, and hard to say out loud. The common-language effect size (CLES) describes the same difference as a probability: pick one participant at random from each condition — how often does the one who could not see outcome variability give the higher estimate? There is a pleasing twist: that quantity has the same probability-of-superiority form as the outcome Zhang et al. asked participants to estimate about the scientific result they saw. The two answer different questions, though: Zhang et al.’s PSup concerns the treatment and control outcomes inside the stimulus; our CLES concerns participants’ PSup judgments across the two visualization conditions.

effectsize’s p_superiority() computes it from the data — the same function the M09 pre-study used. (parametric = FALSE counts every cross-condition pair, ties as half; fct_rev() points it the same direction as your d.)

Before quoting it, check yourself: does a CLES of, say, 80% mean that 80% of the participants in one condition scored higher than every participant in the other?

No. CLES counts pairs, not people, and it only asks who was higher — never by how much. It means: draw one person from each condition, and the one who could not see outcome variability gives the higher estimate about that often. Look back at your section-b figure: the two clouds overlap heavily. A reliable tilt, not a separation — knowing someone’s condition tells you which way to bet, not what number they gave. (To feel what your d means, Dr. Kristoffer Magnusson’s interactive Cohen’s d lets you drag d and watch the overlap and this probability change together. It draws smooth Normal curves, so treat it as intuition, not a picture of your coarse 50–100 responses.)

f · APA result line

Fill this in with your numbers

Among N = ___ [population], those who [saw inferential uncertainty only] estimated a higher probability of superiority (M = ___, SD = ___) than those who [saw outcome variability] (M = ___, SD = ___), a difference of ___ percentage points, 95% CI [___, ___], t(___) = ___, p ___, Cohen’s d = ___, 95% CI [___, ___].

The bracketed phrases come from your study’s row (the paper does not report Cohen’s d for this visualization-condition contrast, so your d is an addition rather than a reproduction target):

Your study [population] [saw inferential uncertainty only] [saw outcome variability] Paper reports
1 · Medical Providers medical providers viewing a hypothetical blood-pressure medication trial first saw the standard-error visualization first saw the standard-deviation visualization t(65.1) = 6.83, p < .001
2 · Data Scientists data scientists were shown standard-error bars only were also shown individual data points t(159.4) = 6.34, p < .001
3 · Faculty tenure-track academic faculty were shown standard-error bars only were also shown individual data points t(363.9) = 4.52, p < .001

Where each number comes from — every one is already on your screen, in the Module’s reporting order: raw effect first, then its uncertainty, then the standardized effect.

Slot Where you got it
N the total across both conditions: nrow(my_data), or the two n values in group_descriptives added together
M, SD per group the M and SD columns of group_descriptives, section b
difference in percentage points estimate from test_result, section d (the printout in c matches it)
95% CI on the difference lower_ci / upper_ci from test_result
tobs, df, p statistic, t_df, p_value from test_result
Cohen’s d and its CI d_result, the cohens_d() output from section e

Two conventions to honor: report p as < .001 rather than as a rounded zero, and keep the direction of your sentence matching the sign of your estimate.

Write it with inline R — the habit Project 2 requires

In your notebook, don’t type these numbers. Let R write them into the sentence — the same two-step pattern you used in the M04 lab and Project 1, now with the APA formatting from your M08 paper (Writing with Code collects the syntax). Then every reported number is linked directly to what your code computed, and if you ever fix something upstream and re-render, the sentence updates itself. Project 2 requires it: anything derived from your data must be generated by code, not typed.

It takes two steps. First, one chunk pulls every number the sentence needs and formats it once. Run it here to check the values, then copy it into a regular {r} chunk in your notebook, right before your APA sentence, and label that chunk apa-numbers:

Three formatting choices sit in that chunk. sprintf() fixes the decimals ("%.2f" is two). p gets the M08 treatment: < .001 when it is that small, and no leading zero otherwise, because a p-value can’t exceed 1. d keeps its leading zero, because it can.

Second, the sentence uses those names with ordinary inline R — expressions like `r apa_N` and `r apa_M_comp`. Paste the template below into your notebook’s APA result section as prose (not inside a chunk), and fill the three bracketed phrases with the words from your study’s row in the table above. Those are words, so you type them; every number comes from R:

Among `r apa_N` [population], those who [saw inferential uncertainty only] estimated a higher probability of superiority (*M* = `r apa_M_comp`, *SD* = `r apa_SD_comp`) than those who [saw outcome variability] (*M* = `r apa_M_ref`, *SD* = `r apa_SD_ref`), a difference of `r apa_diff` percentage points, 95% CI [`r apa_lo`, `r apa_hi`], *t*(`r apa_df`) = `r apa_t`, *p* `r apa_p`, Cohen's *d* = `r apa_d`, 95% CI [`r apa_d_lo`, `r apa_d_hi`].

Inline R is just R placed inside a sentence: `r apa_N` means insert the current value of apa_N here. When your notebook renders, R replaces each expression with the value stored in that object. Check the rendered sentence against your printed outputs — if a number is wrong, fix the code that created it, not the sentence.

Checkpoint 3 · Your study reproduced

For your assigned study you have all six pieces: a data picture showing error bars and points, an estimation plot putting the difference and its uncertainty on their own axis, a Welch’s t-test, a Cohen’s d with its 95% CI, a completed APA result line, and a reproduction verdict in one of three words: exactly reproduced, closely reproduced (with the discrepancy explained), or not reproduced. Your numbers are ready to share.

One check worth making before you present: does your estimation plot agree with your APA line? The difference the plot shows should have the same sign and roughly the same interval as the t-test’s. If the plot says +15 and your sentence says the groups went the other way, one of the two has its conditions reversed.

In your notebook — this whole step is the deliverable core: add the data picture (label the chunk fig-psup with a fig-cap), the estimation plot (label it fig-effect, keep the subtitle key, and write a fig-cap that states the finding in words, without numbers — for example, “[Comparison condition] participants estimated higher PSup than [baseline condition] participants; the bootstrapped difference and its 95% interval sit above zero.” Keep the numbers in your APA sentence: a fig-cap is a chunk option, and chunk options can’t run inline R, so any number there would have to be typed), the Welch’s t-test, the Cohen’s d with its CI, your APA result line — written with inline R, using section f’s apa-numbers chunk and sentence template — and the reproduction verdict. Interpretation prose after each. Your Data section should still say in a sentence or two where the file came from, what one row is, and which rows you analyzed — that is the documentation habit Project 2 will ask you to formalize.

Step 4 · Share back and synthesize

Each group’s spokesperson takes ~2 minutes:

  1. Question — what was the manipulation, in plain English?
  2. Result — your APA line, on screen.
  3. Reproduction verdict — your tobs beside the paper’s reported t. (The paper does not report Cohen’s d for this visualization-condition contrast, so your d is an addition rather than a reproduction target.) Did the numbers come back out?
  4. One sentence — what does your study say about the paper’s central claim?

As each group reports, we’ll fill in this synthesis table on the screen. By the end, we have a 3-row reproduction checklist for the paper’s focal PSup result across all three studies.

Put the three estimation plots side by side

Once all three groups have reported, pull up the three estimation plots together and read the right-hand axes across them. The difference shrinks as you move from Study 1 to Study 3 — a steady decline across the studies that we’ll call the effect-size gradient — and the intervals show how precisely each difference is pinned down.

Two questions worth arguing about before you leave:

  • Do all three effects point the same way? If every interval sits entirely above zero, all three studies agree on the direction: hiding outcome variability raises average PSup judgments.
  • How similar are the magnitudes? Compare the point estimates and interval widths descriptively — and stop there. Whether two intervals happen to overlap is not a test of whether the effects are equal, and non-overlap is not proof that they differ. A formal comparison is a different procedure than eyeballing two pictures.

And there’s a reason this lab cannot settle it anyway. Study 1 differs from Studies 2 and 3 in several ways at once: the audience, the scenario (a medication trial versus a published-study abstract), the visualization contrast (SE vs SD bars, against SE-only vs SE-plus-points), and the design (a first response from a crossover sequence versus straightforward between-subjects assignment). When several things change together, no amount of staring at the gradient can tell you which one produced it — and ordinary sampling variation is tangled in too. Describing the pattern is honest. Explaining it would require a design built for that question.

This is still the payoff of plotting effects rather than p-values: three p-values would all read “< .001” and tell you nothing about how big the illusion is in each audience.

Checkpoint 4 · The focal result, across all three studies

The class synthesis table has all three rows filled — each study’s tobs, d, and verdict beside the paper’s reported values. You can describe the effect-size gradient across the three studies, and explain why this set of studies cannot by itself identify what produced it.

In your notebook — the synthesis table is a class artifact, not required in your submission — but your Discussion draws on it. The Final render section below lists the four things your Discussion needs.

Fill the table in live if you are in class; this is the backstop if you were away or are working ahead. Each row is what that group’s own lab path produces from the shipped data.

Study Course-file N Paper’s reported t Our tobs Our d [95% CI]
1 · Medical Providers 75 6.83 t(65.1) = 6.83 1.60 [1.06, 2.12]
2 · Data Scientists 175 6.34 t(159.2) = 6.41 0.96 [0.64, 1.27]
3 · Faculty 368 4.52 t(363.9) = 4.52 0.47 [0.26, 0.67]

The gradient is in the last column: 1.60 → 0.96 → 0.47. All three reproduce the paper’s direction and significance, and the standardized effect shrinks markedly from Study 1 to Study 3.

Before you write that up as “the illusion is weaker in faculty,” re-read limitation 3 above. The three studies differ in several ways at once — audience, scenario, the visualization contrast, and the design — so a gradient across them cannot say which difference produced it.

Lab debrief · 5 minutes

That’s the reproduction done, and the paper’s focal result assembled across all three studies. With a few minutes left, save your work and look up — this is the last lab of Part 2, so it is worth spending them on what the last four weeks were actually for.

Lab debrief · what did we learn by doing?

  1. The sticking point. Where did today cost you the most time — keeping the sign and order consistent across the plot, the test, and d; reading the estimation plot; or assembling the APA line? What got you past it, and what would have helped sooner?

  2. Running the test was one line — the decisions weren’t. Every group ran the same one-line Welch call. The work sat around it: which rows form the analytic sample (Study 2’s whole discrepancy was one participant), which direction the comparison runs, and how to report it honestly. When you meet a new dataset with no lab page attached, what is the first question you ask?

  3. The same logic, three different audiences. Your group reproduced one study; the other two reproduced the others. The t-values ranged from 6.83 down to 4.52 and the standardized effects from 1.60 down to 0.47 — yet all three would be reported as p < .001. What does that tell you about how much a p-value on its own says about the size of a finding?

  4. What a reproduction is worth — and what it isn’t. You recovered the paper’s numbers from the paper’s data. Say precisely what that establishes and what it leaves open. If a reproduction had failed — your tobs nowhere near their reported t — name three things that could explain it, only one of which is “the paper is wrong.”

  5. Part 2 ends here, and Project 2 starts from it. Everything you did today is the shape of Project 2: find a published result, reconstruct the sample, re-run the test, and report honestly whether it came back. What was hardest today that you would want to have settled before committing to a paper of your own?

Final render and submit

This is the take-home half of the lab. You scaffolded m09_lab.qmd in Step 0 and each step’s In your notebook note told you what to drop in. Now assemble it into a report a stranger could follow cold — and keep building the habit you’ll lean on for Project 2.

What your # Discussion needs — four things

Your Results section had six labelled parts to work through. Discussion has four, and it is the shortest section in the report — aim for one paragraph, not a page.

  1. Did you reproduce it? One sentence, naming your own numbers: “Our reanalysis recovered the published contrast, t(…) = …, d = …” — or saying plainly where it diverged.
  2. Where your study sits. One sentence placing your group against the other two — the effect-size gradient across the three audiences.
  3. One limitation. The strongest candidate is right in front of you: your group is one audience, and the three studies differ in several ways at once — sample, scenario, and what was manipulated — so a gradient across them is not a clean test of who is more susceptible.
  4. What a reproduction does and does not establish. You re-ran their analysis on their data. Agreement shows the published numbers follow from the data as analyzed — computational reproducibility. It is not independent replication, which would need new participants.

The submission steps:

  1. Add your name to the author: field.
  2. Do a final render. Click Render, or press Cmd/Ctrl + Shift + K.
  3. Open the rendered HTML (it renders next to the .qmd in programs/).
  4. Read it end to end — as if a stranger opened it cold. Is every result labeled? Is your APA line complete?
  5. Submit it to Canvas under “Lab 9 — Reproducing the Illusion of Predictability.”
  6. Commit and push. In GitHub Desktop, commit your lab notebook with a one-line summary, then Push origin. Project 2 starts next week and your team will work through this same loop — today is the last low-stakes rehearsal.

One thing to bring to the Project 2 launch. As you assemble the report, note which part of today’s workflow you’d feel ready to scale up to a full reproduction — and which part still feels under-scaffolded. We open next week’s session with that question.

Double check

Before you submit:

What you just did, in research terms

You ran a computational reproduction: you used the data released with a published paper to reconstruct and rerun one of its focal analyses, then compared your result with the published one. (Not a direct replication: that would mean collecting a new sample of clinicians, data scientists, or faculty and running the study again.) That is the credibility check the field increasingly asks for, and it’s the spine of Project 2. Notice what made it trustworthy: you pictured the data before testing it, you chose Welch’s t so the inference did not require equal population variances, you reported an effect size with a confidence interval rather than a lone p-value, and you documented the analytic dataset in your write-up so someone else could retrace your path from raw file to result. Across three expert audiences the illusion held — the same confusion Hofman first found in a boulder-sliding game turns up in clinicians reading a medication-trial scenario. Same recipe, real stakes.

More practice (optional)

After class, any time. Nothing here goes into your lab notebook, and none of it is graded.

  • The same clinician, twice: a paired look at Study 1 → — the paired analysis the lab set aside. Run the paired t-test on both of each clinician’s estimates, then meet two surprises the first-response comparison sidesteps: the within-person correlation is essentially zero, so pairing’s textbook covariance advantage is absent (pairing still buys precision here, just through a different route), and the size of the within-person contrast differs by presentation order. A useful look at what “paired designs are more precise” actually claims.

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