Fix course mismatches found in a full re-audit

src/content/exercises/module-04.ts

@@ -248,7 +248,9 @@ export const module04: ExerciseDef[] = [
               } else if (is.null(labels$x) || is.null(labels$y) ||
                          identical(labels$x, "autonomy") || identical(labels$y, "wellbeing")) {
                 list(pass = FALSE, message = "The axes are still labelled with the column names. A figure has to be readable on its own, by someone who has never seen your data frame: labs(x = \\"Autonomy (1 to 10)\\", y = \\"Wellbeing (0 to 100)\\").")
-              } else if (!inherits(p$theme$panel.grid, "element_blank")) {
+              } else if (!inherits(p$theme$panel.grid, "element_blank") &&
+                         # Older ggplot2 (3.4.4 at least) blanks panel.grid.major and .minor, not panel.grid itself.
+                         !(inherits(p$theme$panel.grid.major, "element_blank") && inherits(p$theme$panel.grid.minor, "element_blank"))) {
                 list(pass = FALSE, message = "The grey panel and its grid lines are still there. APA figures leave out decoration that carries no information, and this exercise asks for the plain white style with axis lines and no grid: add + theme_classic().")
               } else {
                 list(pass = TRUE, message = "That figure could go in a report as it stands: the points show every employee, the line shows the model, the axes say what they mean, and nothing is drawn that carries no information. Module 10 fits exactly this line with lm() and tells you whether its slope is worth believing.")

src/content/exercises/module-07.ts

@@ -186,7 +186,7 @@ export const module07: ExerciseDef[] = [
           value <- as.numeric(value)
           # A band, not a value: a 95 per cent procedure captures mu about 95 times
           # in 100, and the count has an SD of about 2.2, so anything from 86 up is
-          # an honest result. Every wrong answer below lands at 100, 5, 0.95 or ~60.
+          # an honest result. Every wrong answer below lands at 100, 9, 0.91 or 48.
           if (value > 0 && value < 1) {
             list(pass = FALSE, message = paste0("That is the PROPORTION (", round(value, 2), "). The exercise asks for the count out of 100, so use sum() rather than mean()."))
           } else if (isTRUE(all.equal(value, 100, tolerance = 1e-9))) {

src/content/exercises/module-10.ts

@@ -262,6 +262,8 @@ export const module10: ExerciseDef[] = [
       'd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nmodel_wl <- lm(wellbeing ~ workload, data = d)\nbreaks <- quantile(d$workload, probs = c(0, 1/3, 2/3, 1))\nd$third <- cut(d$workload, breaks = breaks, include.lowest = TRUE, labels = FALSE)\nwl_means <- aggregate(wellbeing ~ third, data = d, FUN = mean)',
       // The same table with an extra column of SDs, which is still one row per third.
       'library(dplyr)\nlibrary(broom)\nd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nmodel_wl <- lm(wellbeing ~ workload, data = d)\nwl_means <- d %>%\n  mutate(third = ntile(workload, 3)) %>%\n  group_by(third) %>%\n  summarise(mean_wellbeing = mean(wellbeing), sd_wellbeing = sd(wellbeing), n = n())',
+      // aggregate() with several summaries: the means arrive inside one matrix column.
+      'library(dplyr)\nd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nmodel_wl <- lm(wellbeing ~ workload, data = d)\nd$third <- ntile(d$workload, 3)\nwl_means <- aggregate(wellbeing ~ third, data = d, FUN = function(x) c(mean = mean(x), sd = sd(x)))',
     ],
     check: `
       if (!has_answer("model_wl") || !has_answer("wl_means")) {
@@ -288,12 +290,10 @@ export const module10: ExerciseDef[] = [
         } else if (nrow(wl_means) != 3L) {
           list(pass = FALSE, message = paste0("wl_means has ", nrow(wl_means), " rows. ntile(workload, 3) splits the employees into three groups, so the table should have three."))
         } else {
-          numeric_cols <- vapply(wl_means, is.numeric, logical(1))
-          means_col <- NULL
-          for (nm in names(wl_means)[numeric_cols]) {
-            if (isTRUE(all.equal(as.vector(wl_means[[nm]]), expected, tolerance = 1e-6, check.attributes = FALSE))) means_col <- nm
-          }
-          if (is.null(means_col)) {
+          # numeric_columns(): aggregate() with a FUN returning several summaries
+          # stores them as one matrix column, which a plain column scan misses.
+          found_means <- any(vapply(numeric_columns(wl_means), function(col) isTRUE(all.equal(col, expected, tolerance = 1e-6, check.attributes = FALSE)), logical(1)))
+          if (!found_means) {
             list(pass = FALSE, message = paste0("No column of wl_means holds the mean wellbeing of the three workload thirds, which are ", paste(round(expected, 1), collapse = ", "), ". Check that you split on workload rather than on wellbeing."))
           } else if (b >= 0) {
             list(pass = FALSE, message = "Your model's workload slope is not negative, which contradicts the table: mean wellbeing falls as workload rises. Check which variable you put on the right of the tilde.")

src/content/exercises/module-12.ts

@@ -302,7 +302,7 @@ export const module12: ExerciseDef[] = [
         } else if (isTRUE(all.equal(sort(as.vector(tbl$p.value)), sort(as.vector(none$p.value)), tolerance = 1e-6, check.attributes = FALSE))) {
           list(pass = FALSE, message = "Those p-values are unadjusted. Six comparisons at .05 each give about a 26% chance of at least one false positive if the tests were independent (about 20% here, since they share groups), which is exactly what the adjustment is for. Pass adjust = \\"tukey\\".")
         } else if (!isTRUE(all.equal(sort(as.vector(tbl$p.value)), sort(as.vector(tukey$p.value)), tolerance = 1e-6, check.attributes = FALSE))) {
-          list(pass = FALSE, message = "Your p-values are adjusted, but not with Tukey's method. Tukey is the one designed for all pairwise comparisons after a linear model; Bonferroni is more conservative here. Pass adjust = \\"tukey\\".")
+          list(pass = FALSE, message = "Your p-values are adjusted, but not with Tukey's method. Tukey is the one designed for all pairwise comparisons after a linear model, and it is the one this lesson reports. Pass adjust = \\"tukey\\".")
         } else if (!is.numeric(n_sig) || length(n_sig) != 1L) {
           list(pass = FALSE, message = "n_significant should be a single number: how many rows of pairs_tbl have p.value below .05.")
         } else if (!isTRUE(all.equal(as.numeric(n_sig), as.numeric(exp_n), tolerance = 1e-9, check.attributes = FALSE))) {

src/content/exercises/module-14.ts

@@ -227,6 +227,9 @@ export const module14: ExerciseDef[] = [
       'library(dplyr)\nlibrary(tidyr)\nlibrary(lmerTest)\nd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nlong_d <- d %>%\n  pivot_longer(cols = c(engagement_t1, engagement_t2), names_to = "time", values_to = "engagement") %>%\n  mutate(time = factor(time, levels = c("engagement_t1", "engagement_t2"), labels = c("t1", "t2")))\nm_time <- lmer(engagement ~ time + (1 | employee_id), data = long_d)\nt_lmer <- summary(m_time)$coefficients["(Intercept)", "t value"]\nt_paired <- t.test(d$engagement_t2, d$engagement_t1, paired = TRUE)$statistic',
       // A one-sample test on the time 2 scores, which runs and answers nothing.
       'library(dplyr)\nlibrary(tidyr)\nlibrary(lmerTest)\nd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nlong_d <- d %>%\n  pivot_longer(cols = c(engagement_t1, engagement_t2), names_to = "time", values_to = "engagement") %>%\n  mutate(time = factor(time, levels = c("engagement_t1", "engagement_t2"), labels = c("t1", "t2")))\nm_time <- lmer(engagement ~ time + (1 | employee_id), data = long_d)\nt_lmer <- summary(m_time)$coefficients["timet2", "t value"]\nt_paired <- t.test(d$engagement_t2)$statistic',
+      // The textbook's REML = FALSE: a maximum-likelihood fit, whose t is close to
+      // the paired t but not equal to it.
+      'library(dplyr)\nlibrary(tidyr)\nlibrary(lmerTest)\nd <- read.csv("data/workplace.csv", stringsAsFactors = TRUE)\nlong_d <- d %>%\n  pivot_longer(cols = c(engagement_t1, engagement_t2), names_to = "time", values_to = "engagement") %>%\n  mutate(time = factor(time, levels = c("engagement_t1", "engagement_t2"), labels = c("t1", "t2")))\nm_time <- lmer(engagement ~ time + (1 | employee_id), data = long_d, REML = FALSE)\nt_lmer <- summary(m_time)$coefficients["timet2", "t value"]\nt_paired <- t.test(d$engagement_t2, d$engagement_t1, paired = TRUE)$statistic',
     ],
     alternateSolutions: [
       // The paired test written the other way round: the t flips sign, and the
@@ -251,6 +254,7 @@ export const module14: ExerciseDef[] = [
         reference <- lmerTest::lmer(engagement ~ time + (1 | employee_id), data = long)
         exp_t <- as.vector(coef(summary(reference))[2, "t value"])
         exp_intercept_t <- as.vector(coef(summary(reference))[1, "t value"])
+        ml_t <- as.vector(coef(summary(lmerTest::lmer(engagement ~ time + (1 | employee_id), data = long, REML = FALSE)))[2, "t value"])
         exp_paired <- as.vector(t.test(d$engagement_t2, d$engagement_t1, paired = TRUE)$statistic)
         unpaired <- as.vector(t.test(d$engagement_t2, d$engagement_t1)$statistic)
         if (!inherits(m_time, "merMod")) {
@@ -259,6 +263,8 @@ export const module14: ExerciseDef[] = [
           list(pass = FALSE, message = "Both t statistics should be single numbers.")
         } else if (isTRUE(all.equal(t_lmer, exp_intercept_t, tolerance = 1e-4, check.attributes = FALSE))) {
           list(pass = FALSE, message = "t_lmer is the intercept's t, which tests whether engagement at time 1 differs from zero. The row you want is the one named after the second level of time.")
+        } else if (isTRUE(all.equal(abs(t_lmer), abs(ml_t), tolerance = 1e-4, check.attributes = FALSE))) {
+          list(pass = FALSE, message = paste0("That t (", round(ml_t, 3), ") comes from a maximum-likelihood fit, REML = FALSE. The textbook recommends ML for fixed effects, but the exact match with the paired t-test holds for lmer's default REML fit. Leave REML = FALSE out and t_lmer becomes ", round(exp_t, 3), "."))
         } else if (!isTRUE(all.equal(abs(t_lmer), abs(exp_t), tolerance = 1e-4, check.attributes = FALSE))) {
           list(pass = FALSE, message = paste0("t_lmer is ", round(t_lmer, 4), " but the fixed effect of time has t = ", round(exp_t, 4), "."))
         } else if (isTRUE(all.equal(abs(t_paired), abs(unpaired), tolerance = 1e-6, check.attributes = FALSE))) {

src/content/exercises/module-15.ts

@@ -142,7 +142,7 @@ export const module15: ExerciseDef[] = [
         } else if (!isTRUE(all.equal(b, exp_b, tolerance = 1e-6, check.attributes = FALSE))) {
           list(pass = FALSE, message = paste0("b_wellbeing is ", round(b, 4), " but the wellbeing coefficient is ", round(exp_b, 4), "."))
         } else {
-          list(pass = TRUE, message = paste0("b = ", round(exp_b, 4), " log odds per point of wellbeing. Negative, so higher wellbeing goes with a lower chance of leaving - which is the direction the leaving rates by third showed in the last lesson. Log odds are not readable as they stand; exp() fixes that in Lesson 15-3."))
+          list(pass = TRUE, message = paste0("b = ", round(exp_b, 4), " log odds per point of wellbeing. Negative, so higher wellbeing goes with a lower chance of leaving - which is the direction the leaving rates by third showed in the last lesson. Log odds are not readable as they stand; exp() fixes that, starting with the next exercise."))
         }
       }
     `,

src/content/lessons/08-3-errors-and-power.mdx

@@ -18,7 +18,7 @@ counts it.
 ## Type I: the null really is true
 
 In this population, `programme` is practically unrelated to `exam_score`:
-Psychology and Business average within a tenth of a point of each other. So run
+Psychology and Business average about a tenth of a point apart. So run
 a thousand studies comparing them and count how often the test cries wolf.
 
 <CodeBlock id="typeone" code={`library(dplyr)

src/content/lessons/12-1-two-groups.mdx

@@ -86,7 +86,7 @@ alternative is the Mann-Whitney test, `wilcox.test(wellbeing ~ remote, data = d)
   choices={[
     { text: 'Because t.test() is outdated and no longer recommended', response: 'It is neither. It is a perfectly good function that answers one specific question.' },
     { text: 'Because the same framework then extends to more groups, several predictors, covariates and interactions without learning a new procedure each time', correct: true, response: 'Right. The t-test is the end of its own road; lm() is the beginning of Modules 11 to 15.' },
-    { text: 'Because lm() makes fewer assumptions', response: 'The assumptions are identical - they are the same model.' },
+    { text: 'Because lm() makes fewer assumptions', response: 'It makes the same assumptions as the pooled t-test, because they are the same model. The Welch test that t.test() runs by default actually assumes less: it drops equal variances.' },
     { text: 'Because lm() gives a smaller p-value', response: 'It gives exactly the same p-value, which is the point of this lesson.' },
   ]}
 />

src/content/lessons/15-2-glm-and-log-odds.mdx

@@ -99,6 +99,11 @@ deviance, 514.30 - 451.47 = 62.83, on 2 degrees of freedom.
 `anova(m_left, test = "LRT")` on the model alone answers a different question: it
 tests each predictor in turn, in the order of the formula.
 
+Lesson 15-3 reads them in full, but the next step is one function: `exp()` of a
+log-odds coefficient is an **odds ratio**, the factor by which the odds are
+multiplied for each extra unit of the predictor. Below 1 the odds shrink; above 1
+they grow.
+
 <Exercise id="m15-2-b" />
 
 <Quiz

src/content/lessons/16-3-bayes-factors.mdx

@@ -133,9 +133,10 @@ ttestBF(formula = wellbeing ~ training, data = d)
 lmBF(wellbeing ~ autonomy, data = d)
 ```
 
-Its default prior differs from the one the BIC implies, so its numbers differ
-somewhat. Report which prior you used, whichever tool you use. The free program
-JASP does the same through menus.
+Its default prior differs from the one the BIC implies, and the difference can
+move a result across a label: for remote it gives BF10 = 3.28 (moderate), where
+the BIC gave 1.62 (anecdotal). Report which prior you used, whichever tool you
+use. The free program JASP does the same through menus.
 
 <Quiz
   id="q-bf-meaning"

src/content/lessons/17-1-mediation.mdx

@@ -101,7 +101,8 @@ model <- "
   total := c + a * b
 "
 set.seed(1)
-fit <- sem(model, data = d, se = "bootstrap", bootstrap = 1000)
+# fixed.x = FALSE: without it, some lavaan versions drop a quarter of the draws as failed
+fit <- sem(model, data = d, se = "bootstrap", bootstrap = 1000, fixed.x = FALSE)
 parameterEstimates(fit, boot.ci.type = "perc")
 ```
 

src/content/lessons/17-2-factor-analysis.mdx

@@ -43,7 +43,7 @@ guide.
 
 <Predict
   id="p-rotation"
-  question="Work pressure and manager support are probably related: people with less support may feel more pressure. Should the two factors be allowed to correlate?"
+  question="Work pressure and manager support are probably related: people under more pressure may also get more help from their manager. Should the two factors be allowed to correlate?"
   choices={[
     { text: 'No, factors must always be independent', response: 'That is a choice, called orthogonal rotation (varimax). For psychological constructs it is rarely realistic.' },
     { text: 'Yes, so use an oblique rotation such as promax', correct: true, response: 'Right. If the factors turn out to be uncorrelated, promax gives nearly the same answer as varimax anyway.' },

src/sims/Distribution.test.tsx

@@ -43,6 +43,15 @@ describe('Distribution', () => {
     }
   });
 
+  test('a narrower SD pulls the cut back inside the slider range', () => {
+    render(<Distribution />);
+    fireEvent.change(screen.getByLabelText(/Cut-off/), { target: { value: '100' } });
+    fireEvent.change(screen.getByLabelText(/Standard deviation/), { target: { value: '3' } });
+    // The slider now ends at 70 + 4 * 3 = 82, and the readout follows it.
+    expect((screen.getByLabelText(/Cut-off/) as HTMLInputElement).value).toBe('82');
+    expect(screen.getByTestId('z-score').textContent).toContain('4.00');
+  });
+
   test('the curve is described in words for screen readers', () => {
     render(<Distribution />);
     expect(screen.getByRole('img').getAttribute('aria-label'))

src/sims/Distribution.tsx

@@ -32,7 +32,9 @@ function curvePoints(mu: number, sigma: number): { x: number; value: number }[]
 export default function Distribution() {
   const [mu, setMu] = useState(70);
   const [sigma, setSigma] = useState(10);
-  const [cut, setCut] = useState(70);
+  const [rawCut, setCut] = useState(70);
+  // Kept inside the slider's range, so a narrower SD cannot leave the readout on a cut the slider no longer shows.
+  const cut = Math.min(Math.max(rawCut, mu - 4 * sigma), mu + 4 * sigma);
   const [axis, setAxis] = useState<'raw' | 'z'>('raw');
 
   // One memo, so the curve, the shading and every number below it are read off