Volume 4 · sixteen days

The test is never the first decision.

Design first, model second, p-value last. This volume walks from what a p-value actually means, through the t-test and ANOVA family with planned contrasts, into regression as the general engine — moderation, mediation, and models for binary and ordinal outcomes. You will finish able to name the correct call for a described design, and to report it the way a methods reviewer wants to read it.

days, design to report

bare p-values reported

numbers per effect: estimate, CI, effect size

Day 25 45 minutes · the sentence you must never write

What a p-value is, and the five things it is not

A p-value is the probability of data at least this extreme if the null hypothesis were exactly true and you sampled this way forever. It is a statement about data under an assumption, not about the assumption. Almost every misreading in the literature comes from silently flipping that conditional round.

The misreading	Why it is wrong
`p = .03` means there is a 3% chance the null is true	That is P(H0 \| data). A p-value is P(data \| H0). Flipping a conditional is not a small slip — the two numbers can differ by an order of magnitude.
`p = .06` means there is no effect	It means this sample did not clear an arbitrary threshold. With n = 30 almost nothing clears it. Absence of evidence is not evidence of absence.
`p < .001` means the effect is large	p conflates effect size with sample size. A trivial correlation is p < .001 at n = 5,000.
p tells you the probability of replication	It does not. A study with p = .05 has roughly 50% chance of replicating at the same n.
Two tests, one p < .05 and one p > .05, means the groups differ	Comparing two p-values is not a test of their difference. Test the interaction.

What to report instead

Three numbers per effect, every time: the estimate in the units of your measure (a 3.9-point BDI-II difference), its 95% confidence interval, and a standardised effect size with its own interval (d = 0.62 [0.31, 0.93]). The p-value comes along at the end as a footnote. If your results section reads as a list of p-values, a methods reviewer will ask you to rewrite it.

# The three numbers, in three lines. This is the whole reporting pattern. library(effectsize) t <- t.test(bdi_post ~ group, data = d) # Welch by default es <- cohens_d(bdi_post ~ group, data = d) t$estimate; t$conf.int; es

#> mean in group CBT mean in group control #> 14.2 18.1 #> [1] -6.09 -1.71 #> attr(,"conf.level") [1] 0.95 #> Cohen's d | 95% CI #> ------------------------- #> -0.62 | [-0.94, -0.31]

Do this now · 15 minutes

Find the last three results sentences you wrote (or three from a paper in your field). Rewrite each one as estimate, interval, effect size, then p. Notice which claims stop being defensible once the interval is visible.

Day 26 50 minutes · build the null yourself

Simulate the null and p stops being mysterious

You do not really understand a p-value until you have built one. Shuffle your group labels a few thousand times, compute the difference each time, and you have the distribution of differences that random assignment alone produces. Your observed difference is either ordinary in that distribution or unusual in it. That is all a p-value ever says — and this permutation version needs no distributional assumptions at all.

set.seed(42) obs <- mean(d$bdi[d$group == "CBT"]) - mean(d$bdi[d$group == "control"]) # Break the link between label and outcome 5,000 times null_diffs <- replicate(5000, { shuffled <- sample(d$group) mean(d$bdi[shuffled == "CBT"]) - mean(d$bdi[shuffled == "control"]) }) p_perm <- mean(abs(null_diffs) >= abs(obs)) # two-sided permutation p c(observed = obs, p_permutation = p_perm) hist(null_diffs, breaks = 40, main = "Differences random assignment alone produces") abline(v = c(-obs, obs), lwd = 2, col = "#e26d5c")

#> observed p_permutation #> -3.907 0.0024

Why this is worth 50 minutes

Once you have seen the null distribution as an object you generated, three things become obvious: why larger n narrows it (and so shrinks p for a fixed effect), why p says nothing about effect size, and why one-sided tests must be declared in advance — you can see the tail you would be claiming. Simulation is also the backbone of power analysis in volume 8, so this hour pays twice.

# The parametric version of the same question, for comparison t.test(bdi ~ group, data = d) # And the same logic with bootstrapped CIs instead of a p-value set.seed(42) boot_diffs <- replicate(5000, { s <- d[sample(nrow(d), replace = TRUE), ] mean(s$bdi[s$group == "CBT"]) - mean(s$bdi[s$group == "control"]) }) quantile(boot_diffs, c(.025, .975))

#> Welch Two Sample t-test #> t = -3.12, df = 243.6, p-value = 0.002 #> 95 percent confidence interval: -6.38 -1.44 #> 2.5% 97.5% #> -6.351162 -1.487043

Do this now · 20 minutes

Run the permutation test on your own two-group comparison, then the Welch t-test. Write down both p-values. If they disagree meaningfully, that is information about your distributions — not a reason to pick the smaller one.

Day 27 50 minutes · Welch is the default

The t-test, and the three decisions inside it

R's t.test() defaults to the Welch version, which does not assume equal variances. That default is correct and you should almost never override it: Welch loses essentially nothing when variances are equal and protects you badly-needed error control when they are not. The three decisions that matter are independence, direction, and what you report alongside.

# Independent groups — Welch, the default. Never set var.equal = TRUE by habit. t.test(bdi_post ~ group, data = d) # Paired data: the formula interface will NOT pair for you. t.test(d_wide$bdi_pre, d_wide$bdi_post, paired = TRUE) # Effect sizes with intervals, matching the design library(effectsize) cohens_d(bdi_post ~ group, data = d) # independent hedges_g(bdi_post ~ group, data = d) # small-sample corrected cohens_d(d_wide$bdi_pre, d_wide$bdi_post, paired = TRUE) # within-subject

#> Welch Two Sample t-test #> t = -3.118, df = 243.63, p-value = 0.002034 #> 95 percent confidence interval: -6.38 -1.44 #> #> Cohen's d | 95% CI #> -------------------------- #> -0.62 | [-1.02, -0.22]

Your design	The call
Two independent groups	`t.test(y ~ g, data = d)` — Welch, effect size `cohens_d()`
Same people twice	`t.test(x1, x2, paired = TRUE)` — d with `paired = TRUE`
One group against a fixed value	`t.test(y, mu = 50)` — report the mean and its CI
Two groups, badly skewed, small n	Welch anyway, plus the permutation test from day 26 as a robustness check
Three or more groups	Not a series of t-tests. Day 29.

The mistake reviewers catch

Running a Shapiro-Wilk test, finding p < .05, and switching to Mann-Whitney. Normality tests are hypersensitive at large n and useless at small n — exactly backwards. And a rank test answers a different question (about stochastic dominance, not means), so the switch changes your hypothesis mid-analysis. Judge normality from a residual plot, and choose the test from the design and the quantity you want to claim.

Do this now · 20 minutes

Run your main two-group comparison with t.test() and cohens_d(). Then write the result as a single APA sentence containing the means, the difference, its CI, d with CI, and the p-value last. Keep that sentence — it is your results-section template.

Day 28 50 minutes · residuals, not raw data

Assumptions: which ones matter, and how much

Linear models do not assume your data are normal. They assume the residuals are approximately normal, with roughly constant variance, and that observations are independent. Those three are not equally important: independence violations are fatal, heteroscedasticity is fixable, and mild non-normality at decent n is almost irrelevant thanks to the central limit theorem.

Assumption	How much it matters	What to do
Independence	Fatal. Inflates error rates severalfold.	Fix by design or model it — volume 5.
Linearity (of the modelled relation)	High. A wrong functional form makes coefficients meaningless.	Residual-vs-fitted plot; add a term or transform.
Constant variance	Moderate. Biases standard errors, not estimates.	Welch, robust SEs (`sandwich`), or model the variance.
Normal residuals	Low at n > 30 for inference about means.	Look at the QQ plot; ignore mild curvature.
No influential outliers	Case-dependent, sometimes total.	Cook's distance; report with and without.

m <- lm(bdi_post ~ group + bdi_pre, data = d) library(performance) check_model(m) # the whole diagnostic panel, annotated, in one figure check_normality(m) # of RESIDUALS check_heteroscedasticity(m) check_outliers(m) # Heteroscedasticity present? Keep the model, fix the standard errors. library(lmtest); library(sandwich) coeftest(m, vcov = vcovHC(m, type = "HC3"))

#> OK: residuals appear as normally distributed (p = 0.412). #> Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.017). #> #> t test of coefficients: #> Estimate Std. Error t value Pr(>|t|) #> groupcontrol 3.9182 1.0512 3.7273 0.000239 ***

The habit to build

check_model() on every model you fit, before you look at a single coefficient. It takes four seconds and it catches the two failures that actually invalidate conclusions — wrong functional form and a handful of dominating cases. Screenshot it into your lab notebook so you can show a supervisor what you checked.

Do this now · 20 minutes

Fit your main linear model and run check_model(). Write one sentence per panel: what it shows, and whether you are satisfied. Where you are not, note the remedy rather than deleting data.

Day 29 55 minutes · and why aov() is not enough

Three or more groups: ANOVA, properly specified

An ANOVA asks one deliberately vague question — do these means differ anywhere? — and answers it with an omnibus F. It is a screening device, not a result. The two traps are sums of squares (base aov() gives Type I, which is order-dependent and wrong for unbalanced designs) and stopping at F instead of reporting the comparisons your hypotheses were actually about.

library(afex) # aov_ez wants long data, an id column, and explicit between/within terms. a <- aov_ez(id = "id", dv = "bdi_post", data = d, between = "condition", type = 3) # Type III SS: the default you want a library(effectsize) eta_squared(a, partial = TRUE) # effect size with a CI, always

#> Anova Table (Type 3 tests) #> Response: bdi_post #> Effect df MSE F ges p.value #> 1 condition 2, 243 41.86 8.92 *** .068 <.001 #> #> # Effect Size for ANOVA (Type 3) #> Parameter | Eta2 (partial) | 95% CI #> ----------------------------------------- #> condition | 0.07 | [0.02, 0.13]

Why not just aov()

aov() uses Type I (sequential) sums of squares, so with unequal cell sizes your F for the first factor depends on the order you typed the terms. afex::aov_ez() sets Type III with sum-to-zero contrasts — the combination that makes main effects interpretable in the presence of interactions — and refuses to run if your data are the wrong shape, which is a feature. Use aov() only for perfectly balanced textbook data.

# The omnibus F is a doorway, not a finding. Go straight to the comparisons. library(emmeans) em <- emmeans(a, ~ condition) em # marginal means with CIs pairs(em, adjust = "holm") # all pairwise, multiplicity-corrected eff_size(em, sigma = sigma(a$lm), edf = df.residual(a$lm))

#> condition emmean SE df lower.CL upper.CL #> control 18.1 0.712 243 16.7 19.5 #> CBT 14.2 0.704 243 12.8 15.6 #> Waitlist 17.7 0.718 243 16.3 19.1 #> #> contrast estimate SE df t.ratio p.value #> control - CBT 3.918 0.988 243 3.965 0.0003 #> control - Waitlist 0.402 1.002 243 0.401 0.6888

Do this now · 25 minutes

Run a one-way analysis on your own three-or-more-group variable with aov_ez(), report partial eta-squared with its interval, then produce the marginal means and pairwise contrasts. Write the result without ever saying "the ANOVA was significant" as a standalone sentence.

Day 30 55 minutes · the comparisons you predicted

Planned contrasts beat everything-against-everything

You rarely hypothesise "some means differ". You hypothesise something specific: both treatments beat control, or the effect grows across doses. Contrasts test exactly that — one test, more power, fewer corrections, and a number that maps onto the sentence in your introduction. Weights must sum to zero; that is the only rule.

library(emmeans) em <- emmeans(m, ~ condition) # levels: control, CBT, Waitlist (in this order) contrast(em, list( "treatments vs control" = c(-1, 0.5, 0.5), # weights sum to zero "CBT vs waitlist" = c( 0, 1, -1) ), adjust = "holm") # An ordered predictor: is the change linear across dose? contrast(emmeans(m2, ~ dose), "poly")

#> contrast estimate SE df t.ratio p.value #> treatments vs control -2.160 0.856 243 -2.523 0.0243 #> CBT vs waitlist -3.516 0.995 243 -3.534 0.0010 #> #> contrast estimate SE df t.ratio p.value #> linear 4.812 1.104 196 4.359 <.0001 #> quadratic 0.318 1.061 196 0.300 0.7646

Correction	Use it when
None	One or two preregistered contrasts that constitute your hypotheses.
Holm	Default for a small family of planned tests. Uniformly more powerful than Bonferroni, same guarantee.
Bonferroni	When a reviewer insists, or when the family is tiny and you want transparency.
Tukey	All pairwise comparisons, exploratory. `pairs(em)` uses it by default.
FDR (BH)	Many tests where some false positives are tolerable — screening, not confirmation.

Where people go wrong

Declaring contrasts "planned" after seeing the means. The distinction between planned and post-hoc is about when you decided, and preregistration is what makes the claim checkable — which is why volume 9 day 96 exists. Also: emmeans orders levels as your factor does, so check levels(d$condition) before writing weight vectors, or you will silently test the wrong comparison.

Do this now · 20 minutes

Write your hypotheses as contrast weight vectors before running anything. Check levels(), then run them with adjust = "holm". If you cannot express a hypothesis as weights, it is not yet specific enough to test.

Day 31 55 minutes · interactions are the point

Two factors: the interaction is the hypothesis

In a 2 × 2 you get two main effects and an interaction, and the interaction is usually why you ran the study — it is the claim that an effect depends on something. Two disciplines: plot it before you interpret it, and when it is present, describe the simple effects rather than the main effects, because a main effect averaged over a moderated relationship can describe nobody.

library(afex) a2 <- aov_ez(id = "id", dv = "recall", data = d, between = c("instruction", "delay"), type = 3) a2 eta_squared(a2, partial = TRUE)

#> Anova Table (Type 3 tests) #> Effect df MSE F ges p.value #> 1 instruction 1, 116 12.41 21.08 *** .154 <.001 #> 2 delay 1, 116 12.41 3.62 + .030 .059 #> 3 instruction:delay 1, 116 12.41 9.77 ** .078 .002

library(emmeans); library(ggplot2) em <- emmeans(a2, ~ instruction | delay) # the | means "within each delay" # Simple effects: the instruction effect AT each level of delay contrast(em, "pairwise", adjust = "none") # Is the difference of differences the thing you claimed? Test it directly. contrast(emmeans(a2, ~ instruction * delay), interaction = "pairwise") # The figure, before any of the above afex_plot(a2, x = "delay", trace = "instruction", error = "between") + theme_paper()

#> delay = short: #> contrast estimate SE df t.ratio p.value #> deep - shallow 3.610 0.643 116 5.614 <.0001 #> #> delay = long: #> contrast estimate SE df t.ratio p.value #> deep - shallow 0.472 0.651 116 0.725 0.4699

The classic overclaim

Finding the effect significant in one condition (p = .001) and not in the other (p = .47) and concluding the effect differs between conditions. That comparison is the interaction, and it must be tested as such — two separate tests are not a test of their difference. Here the interaction is significant, so the simple-effects story is licensed. Without it, you report the main effect and stop.

Do this now · 25 minutes

Fit your factorial design, plot the interaction first, then decide from the interaction test whether you are reporting simple effects or main effects. Write the sentence that would appear in the abstract either way.

Day 32 55 minutes · same people, several times

Repeated measures, sphericity, and long format

When each participant contributes several observations, the observations are correlated and the error term has to respect that. Repeated-measures ANOVA does this by partitioning subject variance — at the price of an assumption called sphericity (equal variances of all pairwise differences) and of dropping any participant with a single missing cell. Both limits are why volume 5 exists, but you need this first.

# Long format is mandatory: one row per observation, not per person. long <- d_wide |> pivot_longer(c(t1, t2, t3), names_to = "time", values_to = "score") library(afex) rm <- aov_ez(id = "id", dv = "score", data = long, within = "time") rm # afex applies Greenhouse-Geisser automatically and says so

#> Anova Table (Type 3 tests) #> Response: score #> Effect df MSE F ges p.value #> 1 time 1.72, 84.2 18.31 34.67 *** .258 <.001 #> --- #> Sphericity correction method: GG #> (df are fractional because GG rescales them)

# Mixed design: one between factor, one within factor mix <- aov_ez(id = "id", dv = "score", data = long, between = "condition", within = "time") # Where the change actually is library(emmeans) em <- emmeans(mix, ~ time | condition) contrast(em, "consec", adjust = "holm") # t1→t2, t2→t3 # afex tells you how many participants it dropped. Always read that line. nrow(d_wide); length(unique(mix$data$long$id))

#> Contrasts are set to contr.sum for the following variables: condition #> Warning: Missing values for 7 ID(s), which were removed before analysis: #> 104, 118, 133, 156, 161, 172, 190 #> [1] 253 #> [1] 246

Three things that bite

Wide data. aov_ez() needs long format and an explicit id. Listwise deletion. One missing timepoint removes the whole participant — with attrition that can cost a fifth of your sample, and it is not missing-at-random-safe. Sphericity. With three or more levels, report the corrected test (GG) and say so. A mixed model, which handles unbalanced and missing occasions natively, avoids the second and third entirely.

Do this now · 25 minutes

Reshape your repeated measure to long format, run aov_ez() with within =, and read the dropped-participant warning out loud. Note how many people your design just lost, and keep that number for day 42 when you refit the same data as a mixed model.

Day 33 50 minutes · adjusting for, carefully

ANCOVA: the legitimate use, and the popular misuse

Adding a covariate to a group comparison does two honest things: it reduces error variance (more power) and it adjusts for a baseline difference that randomisation did not eliminate. It does not turn an observational comparison into a causal one, and "controlling for" the wrong variable actively destroys your estimate.

# Randomised pre-post trial: baseline as a covariate is the powerful analysis. m_anc <- lm(bdi_post ~ condition + bdi_pre, data = d) library(car); Anova(m_anc, type = 3) library(emmeans) emmeans(m_anc, ~ condition) # ADJUSTED means, at the mean covariate pairs(emmeans(m_anc, ~ condition), adjust = "holm") # Check the homogeneity-of-slopes assumption explicitly anova(m_anc, lm(bdi_post ~ condition * bdi_pre, data = d))

#> Anova Table (Type III tests) #> Sum Sq Df F value Pr(>F) #> condition 612.4 2 18.312 2.41e-08 *** #> bdi_pre 2841.7 1 169.926 < 2e-16 *** #> #> condition emmean SE df lower.CL upper.CL #> control 17.9 0.512 242 16.9 18.9 #> CBT 14.5 0.508 242 13.5 15.5 #> #> Model 1: bdi_post ~ condition + bdi_pre #> Model 2: bdi_post ~ condition * bdi_pre #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 2 240 3985.1 2 26.41 0.7953 0.4527 ← slopes are homogeneous, fine

Covariate	Verdict
Baseline score on the outcome, randomised trial	Yes. The standard, most powerful analysis of a pre-post RCT.
A nuisance variable measured before treatment	Yes, if preregistered. Otherwise you are choosing covariates by result.
Something caused by the treatment (a mediator)	No. This subtracts the effect you are trying to estimate.
A variable that differs between non-random groups	Dangerous. Lord's paradox: the adjusted comparison may answer no real question.
Anything measured after the outcome	No. It cannot be a cause of the outcome.

Lord's paradox in one sentence

With non-randomised groups, comparing raw change and comparing baseline-adjusted scores can point in opposite directions, and neither is "the" right answer — the correct analysis depends on a causal assumption about how the groups came to differ, which the data cannot supply. Randomisation is what makes the adjusted comparison unambiguous. Without it, draw the causal diagram before choosing the model.

Do this now · 20 minutes

List every covariate you intend to include and write one sentence each justifying it causally — measured before treatment, not affected by it, plausibly related to the outcome. Delete any you cannot justify, and preregister the list.

Day 34 50 minutes · a different question, not a safer one

Rank tests and the bootstrap: what each actually buys

A Wilcoxon test is not a robust t-test. It tests whether one distribution tends to produce larger values than another, which is a different hypothesis with a different effect size, and it cannot adjust for covariates. Often the better robust option is to keep the model you want and get its uncertainty from a bootstrap.

# Rank-based tests, with the effect sizes that belong to them wilcox.test(score ~ group, data = d, conf.int = TRUE) effectsize::rank_biserial(score ~ group, data = d) kruskal.test(score ~ condition, data = d) # 3+ groups effectsize::rank_epsilon_squared(score ~ condition, data = d) # Paired version wilcox.test(d_wide$t1, d_wide$t2, paired = TRUE)

#> Wilcoxon rank sum test with continuity correction #> W = 5824, p-value = 0.0031 #> 95 percent confidence interval: -5.00 -1.00 #> #> r (rank biserial) | 95% CI #> --------------------------------- #> -0.24 | [-0.39, -0.08]

# Keep your model, bootstrap its uncertainty. Works for anything. library(boot) med_diff <- function(data, i) { s <- data[i, ] median(s$score[s$group == "a"]) - median(s$score[s$group == "b"]) } set.seed(1); b <- boot(d, med_diff, R = 5000) boot.ci(b, type = "bca") # Bootstrapped CIs for any lm coefficient, no normality assumption needed library(car) Boot(lm(score ~ group + age, data = d), R = 2000) |> confint(type = "bca")

#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS #> Level BCa #> 95% (-5.00, -1.00 ) #> #> 2.5 % 97.5 % #> (Intercept) 41.12083 47.33914 #> groupb -4.61207 -1.02884 #> age -0.07219 0.11842

Choosing between them

Ask what you want to claim. If the claim is about means — because your theory is about magnitudes and you need covariates — keep the linear model and bootstrap it. If the claim is genuinely about which group tends to score higher, with ordinal data and no covariates, the rank test is the right tool and rank-biserial r is its effect size. Never switch tests because of a normality p-value.

Do this now · 20 minutes

Take one comparison and run three analyses: Welch t-test, Wilcoxon with rank-biserial r, and a bootstrapped mean difference. Write one sentence on what each estimates. Choose the one that matches your hypothesis, and mention the others as robustness checks.

Day 35 45 minutes · the workhorse, and its traps

Correlation: intervals, matrices and attenuation

Correlations are easy to compute and easy to overread. Three habits fix most of it: always report the confidence interval (an r of .30 at n = 40 has an interval from .00 to .55 — a range covering "nothing" to "substantial"), decide Pearson versus Spearman from the shape of the data, and remember that unreliable measures shrink correlations toward zero.

cor.test(d$agree, d$neuro) # r with CI and p cor.test(d$agree, d$neuro, method = "spearman") # monotone, rank-based library(psych) corr.test(d[, c("neuro", "agree", "extra", "consc")], adjust = "holm") # matrix with n and adjusted p print(corr.test(d[, 1:4])$ci, digits = 2) # the intervals, which matter most partial.r(d[, c("neuro", "agree", "age")], c(1, 2), 3) # partial r

#> Pearson's product-moment correlation #> t = -12.94, df = 2434, p-value < 2.2e-16 #> 95 percent confidence interval: -0.288 -0.214 #> sample estimates: cor -0.2513 #> #> raw.lower raw.r raw.upper #> neuro-agree -0.29 -0.25 -0.21 #> neuro-extra -0.24 -0.20 -0.16

Attenuation, and why it matters for your thesis

If two scales have reliabilities of .70, the maximum correlation you can observe between the constructs is about .70 — not 1.00 — because measurement error is not shared. An r of .35 between two such scales corresponds to a disattenuated correlation near .50. This is not licence to report the corrected value as your result; it is the reason to report reliabilities alongside every correlation matrix, and the reason volume 6 comes before you interpret any of this.

# Visualise before interpreting: r is only meaningful for a linear pattern library(GGally) ggpairs(d[, c("neuro", "agree", "extra", "consc")]) # A publishable correlation table straight to Word (volume 9 day 87) apaTables::apa.cor.table(d[, c("neuro","agree","extra")], filename = "output/tables/table1.doc")

Do this now · 15 minutes

Build your correlation matrix with corr.test(), then print the interval matrix. Find the correlation whose interval crosses zero widest and write down the n you would need to pin it to ±.10. That number is your day-78 power analysis, arriving early.

Day 36 55 minutes · the engine underneath everything

lm(): one function that contains most of this volume

The t-test, ANOVA and ANCOVA are all special cases of the linear model. Once you read lm() output fluently you have a single framework for continuous and categorical predictors together, which is how real behavioural questions arrive. The skill is reading each coefficient as a sentence.

d <- d |> mutate(condition = fct_relevel(condition, "control")) # reference level m <- lm(bdi_post ~ condition + bdi_pre + age, data = d) summary(m) confint(m) # intervals for every coefficient effectsize::standardize_parameters(m) # betas, done correctly

#> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 6.412034 1.884120 3.403 0.000771 *** #> conditionCBT -3.407881 0.716402 -4.756 3.39e-06 *** #> conditionWaitlist -0.318442 0.724119 -0.440 0.660450 #> bdi_pre 0.681204 0.052260 13.035 < 2e-16 *** #> age -0.041118 0.031004 -1.326 0.185999 #> #> Residual standard error: 4.07 on 242 degrees of freedom #> Multiple R-squared: 0.526, Adjusted R-squared: 0.518 #> F-statistic: 67.1 on 4 and 242 DF, p-value: < 2.2e-16

Output	Read it as
`(Intercept)`	Predicted outcome when every predictor is zero and each factor is at its reference level. Meaningless unless you centred — which is why day 38 centres.
`conditionCBT -3.41`	CBT scores 3.41 points lower than control, holding baseline and age constant.
`bdi_pre 0.68`	One point higher at baseline predicts 0.68 points higher after, within condition.
`Std. Error`	The precision of that estimate. Estimate ± 1.96 SE is roughly the 95% CI.
`R-squared 0.526`	Share of outcome variance the model accounts for in this sample. It always rises when you add predictors — see day 99.

Dummy coding, said plainly

A k-level factor becomes k − 1 indicator columns, and every coefficient is a comparison against the reference level. So fct_relevel() is a substantive decision, not tidying: put your control or comparison group first and the coefficients become the contrasts you wanted to report. If the intercept is uninterpretable, you have found the reason people centre continuous predictors.

Do this now · 25 minutes

Fit your main model with lm(), set the reference level deliberately, and write one plain sentence per coefficient including its interval. Then produce standardised betas and note which coefficients change rank order — that tells you how different your predictor scales are.

Day 37 50 minutes · before you believe a coefficient

Diagnostics: influence, collinearity, and misspecification

A model can fit well, print three asterisks, and be driven by four participants. The diagnostics worth running are about influence (which cases determine the estimate), collinearity (whether your predictors can be separated at all) and functional form (whether the straight line was right).

library(performance); library(car) check_model(m) vif(m) # >5 worth a look, >10 a real problem check_collinearity(m) infl <- influence.measures(m) which(apply(infl$is.inf, 1, any)) # flagged cases plot(m, which = 5) # residuals vs leverage, Cook contours

#> # Check for Multicollinearity #> Low Correlation #> Term VIF VIF 95% CI Tolerance #> condition 1.04 [1.01, 1.31] 0.96 #> bdi_pre 1.12 [1.05, 1.36] 0.89 #> age 1.09 [1.03, 1.32] 0.92 #> [1] 17 88 204

# The honest response to influential cases: report both fits. m_all <- lm(bdi_post ~ condition + bdi_pre, data = d) m_excl <- lm(bdi_post ~ condition + bdi_pre, data = d[-c(17, 88, 204), ]) library(modelsummary) modelsummary(list("All cases" = m_all, "Excluding 3 flagged" = m_excl), statistic = "conf.int", stars = FALSE)

#> +------------------+-------------------+----------------------+ #> | | All cases | Excluding 3 flagged | #> +------------------+-------------------+----------------------+ #> | conditionCBT | -3.408 | -3.201 | #> | | [-4.819, -1.997] | [-4.552, -1.850] | #> | bdi_pre | 0.681 | 0.694 | #> | | [ 0.578, 0.784] | [ 0.593, 0.795] | #> | Num.Obs. | 247 | 244 | #> +------------------+-------------------+----------------------+

The rule about outliers

Delete a case only for a documented reason that is independent of its effect on your result — a failed attention check, an impossible value, a known equipment fault — and state the rule in your preregistration. Otherwise keep it and show a sensitivity analysis like the table above. "The effect holds with and without the three most influential cases" is a strong sentence; a quietly trimmed dataset is a finding nobody can trust.

Do this now · 20 minutes

Run check_model() and vif() on your model, identify any flagged cases, and build the two-column sensitivity table. Decide your exclusion rule in writing before looking at which direction the estimate moves.

Day 38 55 minutes · does the effect depend on X

Moderation: product terms, centering, simple slopes

Moderation is an interaction between continuous predictors, and it is where interpretation goes wrong most often. Three rules: centre (or use meaningful zero points) so the main effects are interpretable, never interpret a product term from its coefficient alone, and always plot it. The coefficient tells you how the slope changes per unit of the moderator — a quantity almost nobody reads correctly from a table.

d <- d |> mutate(stress_c = as.numeric(scale(stress, scale = FALSE)), support_c = as.numeric(scale(support, scale = FALSE))) m_mod <- lm(depression ~ stress_c * support_c + age, data = d) summary(m_mod) effectsize::eta_squared(car::Anova(m_mod, type = 3), partial = TRUE)

#> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 14.82104 0.31182 47.531 < 2e-16 *** #> stress_c 0.48213 0.06104 7.899 3.12e-13 *** #> support_c -0.31904 0.07211 -4.425 1.36e-05 *** #> stress_c:support_c -0.09817 0.02988 -3.285 0.00114 ** #> #> ← interpretation: each SD more support reduces the stress slope by ~0.10

library(marginaleffects); library(ggeffects) # Simple slopes at conventional moderator values slopes(m_mod, variables = "stress_c", newdata = datagrid(support_c = c(-1, 0, 1) * sd(d$support_c))) # Johnson-Neyman: where does the stress effect stop being distinguishable from 0? interactions::johnson_neyman(m_mod, pred = stress_c, modx = support_c) # The figure. Always. ggpredict(m_mod, terms = c("stress_c", "support_c [-1.5, 0, 1.5]")) |> plot()

#> support_c Estimate Std. Error z Pr(>|z|) 2.5 % 97.5 % #> -2.911 0.768 0.091 8.44 <0.001 0.590 0.946 #> 0.000 0.482 0.061 7.90 <0.001 0.362 0.602 #> 2.911 0.196 0.084 2.33 0.020 0.031 0.361 #> #> JOHNSON-NEYMAN INTERVAL #> The slope of stress_c is significant outside the interval [3.91, 7.44].

Two traps

Uncentred predictors. Without centering, the "main effect" of stress is its slope when support equals zero — often outside the observed range, occasionally impossible. The interaction coefficient itself is unaffected, but everything else becomes uninterpretable. Power. Detecting an interaction typically needs several times the sample size of the corresponding main effect; a null interaction at n = 80 is uninformative, not evidence of no moderation.

Do this now · 25 minutes

Centre your predictors, fit the product term, and produce the simple-slopes figure with ggpredict(). Write the interaction as a sentence containing a number: "each additional point of support reduced the stress-depression slope by 0.10 [0.04, 0.16]".

Day 39 55 minutes · and its causal small print

Mediation: bootstrapped indirect effects, honestly framed

Mediation asks whether X affects Y through M. The statistics are straightforward — a product of two coefficients with a bootstrapped interval — and the causal assumptions are severe: no unmeasured confounding of either path, correct temporal order, and no measurement error in M. Cross-sectional mediation satisfies none of these, which is why reviewers have become sceptical of it.

library(lavaan) model <- ' # paths M ~ a*X Y ~ b*M + cprime*X # derived quantities indirect := a*b total := cprime + (a*b) prop_med := (a*b) / (cprime + a*b) ' fit <- sem(model, data = d, se = "bootstrap", bootstrap = 5000) parameterEstimates(fit, boot.ci.type = "bca.simple", standardized = TRUE)

#> lhs op rhs label est se z pvalue ci.lower ci.upper std.all #> 1 M ~ X a 0.412 0.061 6.754 0.000 0.294 0.533 0.389 #> 2 Y ~ M b 0.318 0.058 5.483 0.000 0.205 0.433 0.301 #> 3 Y ~ X cprime 0.104 0.062 1.677 0.094 -0.017 0.226 0.098 #> 7 ab := a*b indirect 0.131 0.031 4.226 0.000 0.077 0.198 0.117

Sobel is obsolete; so is the Baron and Kenny sequence

The product ab is not normally distributed, so the Sobel test is mis-calibrated — use a bootstrapped BCa interval, as above. And a significant total effect is not a prerequisite for mediation: suppression and competing paths can produce a real indirect effect with no total effect. Report ab with its interval as the result; c' is not "the remaining direct effect" unless the causal assumptions hold.

# Multiple mediators, in parallel, with contrasts between them model2 <- ' M1 ~ a1*X M2 ~ a2*X Y ~ b1*M1 + b2*M2 + cprime*X ind1 := a1*b1 ind2 := a2*b2 diff := ind1 - ind2 # which mediator carries more of the effect ' fit2 <- sem(model2, data = d, se = "bootstrap", bootstrap = 5000) # The causal-inference framing, with sensitivity to unmeasured confounding library(mediation) med <- mediate(lm(M ~ X + age, d), lm(Y ~ X + M + age, d), treat = "X", mediator = "M", boot = TRUE, sims = 2000) summary(med); summary(medsens(med)) # how strong must a confounder be to kill it

#> Estimate 95% CI Lower 95% CI Upper p-value #> ACME 0.131 0.077 0.198 <2e-16 #> ADE 0.104 -0.017 0.226 0.094 #> Total Effect 0.235 0.108 0.361 <2e-16 #> Prop. Mediated 0.557 0.312 0.918 <2e-16

Do this now · 25 minutes

Fit your mediation in lavaan with bootstrapped BCa intervals. Then write the limitations sentence first: name the temporal design, the confounders you did not measure, and the reliability of your mediator. If that sentence undermines the analysis, say so in the paper rather than hoping.

Day 40 55 minutes · when the outcome is not continuous

Binary and ordinal outcomes: glm, odds ratios, and Likert items

Relapse or not, correct or incorrect, five response categories on a single item — none of these are continuous, and a linear model on them produces impossible predictions and wrong standard errors. glm(family = binomial) handles binary outcomes; cumulative link models handle ordered categories. The reporting skill is converting coefficients into something a reader can picture.

m_log <- glm(relapse ~ condition + severity + age, data = d, family = binomial) summary(m_log) exp(cbind(OR = coef(m_log), confint(m_log))) # odds ratios with CIs # Odds ratios are hard to intuit. Predicted probabilities are not. library(marginaleffects) avg_comparisons(m_log, variables = "condition") # risk difference, in probability library(ggeffects); ggpredict(m_log, terms = "severity [all]") |> plot()

#> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> conditionCBT -0.8412 0.2914 -2.887 0.00389 ** #> severity 0.3106 0.0782 3.972 7.13e-05 *** #> #> OR 2.5 % 97.5 % #> conditionCBT 0.4313 0.2417 0.7621 #> severity 1.3643 1.1713 1.5932 #> #> Term Contrast Estimate Std. Error z 2.5 % 97.5 % #> condition CBT - control -0.163 0.055 -2.95 -0.271 -0.055 #> ← 16 percentage points lower relapse probability

# A single Likert item is ordinal, not interval. Model it as such. library(ordinal) d$agreement <- factor(d$agreement, ordered = TRUE) # 1..5 m_ord <- clm(agreement ~ condition + age, data = d) summary(m_ord) nominal_test(m_ord) # is the proportional-odds assumption tenable? # Counts (number of sessions attended, errors made) m_cnt <- glm(sessions ~ condition, data = d, family = poisson) performance::check_overdispersion(m_cnt) # if overdispersed: family = quasipoisson # or MASS::glm.nb()

#> Coefficients: #> Estimate Std. Error z value Pr(>|z|) #> conditionCBT 0.7214 0.2183 3.305 0.00095 *** #> Threshold coefficients: #> Estimate Std. Error z value #> 1|2 -2.1043 0.2311 -9.105 #> 2|3 -0.5128 0.1802 -2.846 #> #> # Overdispersion test #> dispersion ratio = 2.418 #> Pearson's Chi-Squared = 596.7 p-value = < 0.001 ← use quasipoisson or NB

The everyday version of this mistake

Averaging a five-point Likert item across participants and running a t-test. For a scale of several items the interval approximation is usually defensible (and volume 6 tells you when); for a single item with a floor, a ceiling and unequal category spacing it is not. Use clm(), report the odds of endorsing a higher category, and you will have pre-empted a reviewer's first comment.

Do this now · 25 minutes

Fit a logistic model on any binary outcome you have, report both the odds ratio and the average risk difference in percentage points, and plot the predicted probabilities. If you have a single-item ordinal outcome, refit it with clm() and compare the conclusions.

Reference

From described design to the correct call

Start from the sentence describing your design, never from the test name you remember.

Design	The call
Two independent groups, continuous outcome	`t.test(y ~ g)` + `cohens_d()`
Same people twice	`t.test(x1, x2, paired = TRUE)` or `lm` on change scores
3+ independent groups	`aov_ez(between =)` then `emmeans` contrasts
3+ occasions, same people	`aov_ez(within =)`, or `lmer` (volume 5) if unbalanced
Groups × occasions	`aov_ez(between =, within =)`; interaction is the hypothesis
Randomised pre-post trial	`lm(post ~ condition + pre)`, then `emmeans`
Two continuous variables	`cor.test()`, with the interval reported
Several predictors, continuous outcome	`lm()`, with `check_model()`
Effect depends on a continuous variable	`lm(y ~ x * mod)`, centred; simple slopes
Effect operates through a variable	`lavaan::sem()` with bootstrapped ab
Binary outcome	`glm(family = binomial)`; report risk difference too
Single Likert item as outcome	`ordinal::clm()`
Counts	`glm(family = poisson)`, check overdispersion
Clustered or repeated observations	`lmer()` — volume 5

Reference

Sentences you can paste and adapt

t-test

CBT (M = 14.2, SD = 5.1) scored lower than control (M = 18.1, SD = 5.4), a difference of 3.9 points, 95% CI [1.7, 6.1], d = 0.62 [0.31, 0.94], t(243.6) = 3.12, p = .002.

ANOVA + contrast

Condition affected post-treatment scores, F(2, 243) = 8.92, p < .001, partial eta-squared = .07 [.02, .13]. The planned treatments-versus-control contrast was 2.16 points, 95% CI [0.48, 3.84], p = .024.

Regression

Controlling for baseline and age, CBT predicted 3.41 points lower symptoms, 95% CI [1.99, 4.82], beta = -0.31, t(242) = 4.76, p < .001. The model accounted for 52.6% of variance.

Moderation

Support moderated the stress-depression association, b = -0.098, 95% CI [-0.157, -0.040], p = .001: the stress slope fell from 0.77 [0.59, 0.95] at low support to 0.20 [0.03, 0.36] at high support.

Mediation

The indirect effect of X on Y through M was 0.131, bootstrapped 95% BCa CI [0.077, 0.198] from 5,000 resamples. Given the cross-sectional design this is consistent with, but not evidence of, a causal pathway.

Logistic

CBT was associated with lower odds of relapse, OR = 0.43, 95% CI [0.24, 0.75], p = .004, corresponding to 16.3 percentage points lower predicted probability, 95% CI [5.5, 27.1].

Reference

Install these once

install.packages(c( "afex", # ANOVA done right: Type III, GG correction, clear output "emmeans", # marginal means and contrasts for any model "effectsize", # effect sizes with confidence intervals "performance", # check_model() and friends "car", # Anova(), vif(), Boot() "marginaleffects", # simple slopes, risk differences, predictions "ggeffects", # model-implied values ready to plot "lavaan", # mediation and, in volume 6, SEM "ordinal", # cumulative link models for Likert outcomes "boot", # general-purpose bootstrap "modelsummary", # side-by-side model tables "interactions" # Johnson-Neyman ))

Checkpoint

Eight questions before volume 5

Next: volume 5

Everything in this volume assumed one independent observation per row. Diary entries, repeated trials, pupils in classrooms and patients within therapists break that assumption — and mixed-effects models fix it. Twelve days from the ICC to growth curves and a full ESM analysis.

Start day 41 →