Multiple Linear Regression

Pre-Study · Module 11 · Fri Nov 6

Welcome to the M11 pre-study. You’ve already read the M11 Module, which extends M10’s single-predictor line to multiple regression: residualized gains, adjusted comparisons, centering, and the precision gain from adjusting for baseline benevolence. You’ve also seen adjusted means with their confidence and prediction intervals, site added as a predictor, the nested-model F-test, and the idea of starting with an estimand — the quantity that would answer the research question. This pre-study gives you guided practice with those ideas. This is where it comes together before Monday.

You’ll watch three short videos and work through guided practice activities after each one. Every activity uses analysis, the 3,965 REACH participants (of 4,598) with complete data on the primary model variables. It’s already loaded in your R sandbox.

As in M10’s pre-study, every activity stays on benevolence. The conceptual steps here are denser — residualized gain, partial effects, confidence versus prediction intervals, the Table 2 fallacy — so we keep the data fixed and walk through the new ideas carefully.

How this page is organized

Three short videos reinforce the Module’s hardest conceptual steps. After each video, practice activities let you do the analysis yourself. Each video’s section closes with a Quick Check — a short set of questions about the concepts you just practiced.

Every code activity uses a three-tab panel:

  • ✍️ Your Code — the starter with blanks for you to fill in
  • 💡 Hint — a nudge in the right direction
  • 👀 Spoiler — the full working code if you get stuck

Work through the ✍️ tab first. Only open 💡 or 👀 after you’ve tried it yourself — struggling a little is how this stuff sticks. The main data are already loaded in a tibble called analysis (n = 3,965), with centered baseline benevolence t1_benevolence.c pre-computed.


Video 1 — From M10 to MLR: Residualized Gain

Video looks blurry? Click the gear icon at the lower right of the player and set Quality to 1080p, or the highest option offered.

What to listen for:

  • A regression’s residual is the observed outcome minus the model’s prediction; in this baseline-to-follow-up design it’s a residualized gain — the part of T2 that baseline benevolence couldn’t account for
  • When treatment-arm residuals from a baseline-only model sit systematically above the control-arm residuals, that gap is — to a very close approximation — the treatment effect adjusted for baseline benevolence
  • The two-step recipe (fit, residualize, group, compare) and the one-step MLR (lm(t2_benevolence ~ t1_benevolence + treat)) land on the same number to three decimal places
  • The MLR shortcut: one regression call automatically does what we’d otherwise have to walk through by hand
  • Reading the MLR coefficients as comparisons between people: comparing two participants with the same baseline benevolence, one in each arm, the model predicts the treatment participant scores about 0.503 points higher at T2. Because the workbook was randomly assigned, this contrast has a causal basis; because our classroom analysis includes only participants with observed T2 scores, that causal reading also depends on assumptions about the missing follow-ups. The baseline slope, 0.585, compares two people in the same arm one point apart at baseline, not a change within one person

Activity 1.1 — Fit the baseline-only model

Start where M10 left off. Fit a simple linear regression of follow-up benevolence on baseline benevolence — no treatment, just T1 predicting T2 — on the whole complete-case analysis sample: both arms and all six sites, not M10’s Hong Kong control arm.

Your task: fill two blanks — the outcome (left of ~) and the predictor (right of ~).

The outcome is the 2-week follow-up score; the predictor is the baseline benevolence score. The variable prefix is t1_ for baseline (T1) and t2_ for follow-up (T2).

You should see two rows — (Intercept) and t1_benevolence. The slope is positive and precisely estimated: people one point higher in baseline benevolence are predicted to be about three-fifths of a point higher at T2. That line is the baseline-only expectation the rest of this section builds on.

Activity 1.2 — Residuals by treatment

Now look at what the baseline-only model can’t predict. Use augment() to attach each participant’s residual to the data, then group by treatment and compare the mean residuals.

Your task: fill three blanks — the model object to augment, the column to group by, and the residual column to summarize.

The model object is the one you fit in Activity 1.1. The grouping variable is the treatment indicator. The column augment() returns for the residual is .resid.

The control-arm average residual is negative and the treatment-arm average residual is positive — nearly mirror images of each other, each about a quarter of a point from zero. The difference of those two means is the residualized treatment contrast — about 0.503 of a benevolence point. On average, treatment-arm participants exceed their baseline-only expectation, while control-arm participants fall a little short of theirs.

Activity 1.3 — Recover the same number with MLR

The two-step procedure (fit, residualize, compare) is conceptually clean, but it has two quirks. First, it residualizes only the outcome, not treatment, so it only approximates the MLR coefficient (random assignment makes treatment nearly unrelated to baseline benevolence, which is why the approximation is so close here). Second, if we regressed those residuals on treat, that second regression wouldn’t know the first one already used a degree of freedom. The fix is to do both steps in one regression, with both predictors on the right-hand side of the formula. Fit the MLR and look at the coefficient on treat.

Your task: fill three blanks — the outcome, the baseline-benevolence predictor, and the treatment predictor.

The two predictors go on the right side of the formula, separated by a +; use the raw baseline benevolence score, t1_benevolence (Activity 2.1 swaps in the centered one). The order doesn’t matter for the fit; conventionally we list the covariate first and the focal predictor second.

Look at the row for treatTreatment. The coefficient is about 0.503 — matching the residual gap from Activity 1.2 to three decimal places. The MLR gives the exact baseline-adjusted treatment contrast in one pass, with the right degrees of freedom built in. The outcome-only residual gap from Activity 1.2 was an extremely close preview because randomization made treatment nearly unrelated to baseline benevolence.

Quick Check

Answer each question — you’ll see green (correct) or pink (incorrect) feedback as you type.

1. A regression’s residual is best described as:

2. True or False: In an RCT, because random assignment makes treatment nearly unrelated to baseline benevolence, the difference of mean residuals from a baseline-only regression should be roughly the same as the MLR coefficient on treatment.

3. Suppose you regressed the residuals from Activity 1.2 on the treatment indicator, something like lm(.resid ~ treat). Compared to lm(t2_benevolence ~ t1_benevolence + treat, data = analysis), what’s the issue with this two-step approach?


Video 2 — Reading the Baseline-and-Treatment Model: Partial Effects, Centering, and the Precision Gain

Video looks blurry? Click the gear icon at the lower right of the player and set Quality to 1080p, or the highest option offered.

What to listen for:

  • Centering a continuous predictor (subtracting its mean) doesn’t change this additive model — the fitted values, slopes, R², and \(\sigma\) stay the same; the intercept’s value and its meaning change
  • The partial-effect rule for the additive linear models in M11: each slope is the predicted difference in \(\hat{y}\) between people one unit apart on that predictor but with the same values on all other predictors. Once a model includes interactions or nonlinear terms, some coefficients need more careful handling.
  • Inference carries over from M10: for each slope, \(H_0\): \(\beta_j = 0\) and \(t_\text{obs} = b_j / SE(b_j)\), judged against one \(t^*\) from the residual df (\(n - k - 1\) = 3,962, with k = 2 slopes); the overall F-test (critical value \(F^*\) ≈ 3.00) asks whether both slopes are 0 at once, and glance() reports it beside R² and \(\sigma\)
  • Precision gain in this dataset: adjusting for a strong baseline covariate barely moves the treatment estimate and substantially reduces its standard error — and regulatory guidance on trial analysis (EMA, FDA) recommends this for exactly this reason

Activity 2.1 — Fit the centered baseline-and-treatment model

Refit the model with the centered baseline benevolence (t1_benevolence.c) in place of the raw baseline benevolence score. The slopes stay the same; only the intercept moves.

The first lines of the code are filled in for you, and they show how the centered column is built: mutate() subtracts the sample mean of t1_benevolence from every participant’s baseline benevolence score. The summarize() check confirms it worked. The centered column averages 0 (R may print a tiny number such as 1e-16, which is 0 up to computer rounding), and its SD matches the raw column’s, because centering moves zero without changing the spread.

Your task: fill one blank — the centered baseline-benevolence column.

The centered predictor’s name ends in .c — it’s the column the mutate() line at the top of the code just built.

The average baseline benevolence score is about 3.15, so a participant who started there has a centered score of 0. The intercept is now about 3.207 — the predicted T2 for the reference-profile participant (Control arm, at the average baseline benevolence). The slope on baseline benevolence and the treat slope are exactly the values from Activity 1.3. In this additive model, centering didn’t change the fitted values or the slopes; it just relocated the intercept to a meaningful in-data anchor.

Activity 2.2 — Inference for each slope, and for the model as a whole

Video 2 laid out the ingredients of regression inference: for each slope, \(H_0\): \(\beta_j = 0\) and \(t_\text{obs} = b_j / SE(b_j)\); one critical value, \(t^*\), from the residual degrees of freedom (\(n - k - 1\), where k counts the slopes, the non-intercept coefficients); and, for the model as a whole, the overall F-test with its own critical value, \(F^*\). This is the inference you learned in M10, extended to two slopes. Rebuild each piece yourself.

Your task: fill three blanks — the degrees of freedom for \(t^*\), the numerator degrees of freedom for \(F^*\) (one per slope), and the denominator of the test statistic.

  • Blank 1: m1_df — the residual degrees of freedom you pulled on the line above.
  • Blank 2: 2 — the overall F-test has one numerator degree of freedom per slope, and this model has two.
  • Blank 3: std.error — every t statistic in a regression is the estimate divided by its standard error.

The ingredients. df = 3,965 − 2 − 1 = 3,962, so \(t^*\) = 1.961 for every slope and \(F^*\) = 2.998 for the overall test. The t-tests split \(\alpha = .05\) across two tails, which is why qt() uses p = 0.975; the F-test uses only the upper tail, so qf() uses p = 0.95.

Each slope. Every by-hand column matches its tidy() neighbor exactly. For the treatment row, \(t_\text{obs}\) = 0.5031 ÷ 0.02608 = 19.29, far beyond \(t^*\), and the 95% CI is 0.503 ± 1.961 × 0.0261 = [0.452, 0.554]. Among participants with the same baseline benevolence, those assigned to the workbook are estimated to score about half a point higher at T2. The slope on baseline benevolence, 0.585 with a 95% CI of [0.560, 0.609], is large and precisely estimated. (When tidy() prints a p-value of 0 because an extremely small value has rounded or underflowed to zero, report it as p < .001. The intercept row tests whether the reference profile’s predicted T2 is 0, which isn’t a question this study asks.)

The model as a whole. glance() reports R² = .389, so the model accounts for about 38.9% of the variation in T2, and \(\sigma\) = 0.820, the typical size of a residual. Its \(F_\text{obs}\) = 1260.9, on 2 and 3,962 degrees of freedom, is far above \(F^*\): we can reject the null hypothesis that both slopes are 0. The F-test says the slopes aren’t all 0; R² says how much variation they account for together.

Activity 2.3 — The precision gain

Now demonstrate the precision-gain story directly. Fit a treat-only model on the same analysis rows (no adjustment for baseline benevolence) and compare its standard error on treat to the adjusted model’s.

Your task: fill two blanks — the outcome and the treatment predictor in the treat-only formula.

A “treat-only” model is the simplest possible RCT regression — outcome on the left, the single treatment indicator on the right, no other predictors.

The point estimate barely moves — 0.519 in the treat-only model against 0.503 in the adjusted model, a shift of less than 0.02. But the adjusted model’s SE on treat is roughly 19% smaller than the treat-only model’s. The sigma column shows why: the leftover noise shrinks, from \(\sigma\) = 1.02 in the treat-only model to 0.82 once baseline benevolence is in the model, and the SE on treat shrinks with it. The statistic column shows the payoff for the test: with nearly the same estimate and a smaller SE, \(t_\text{obs}\) grows from 16.07 to 19.29. That’s the precision gain. Under randomization, adjusting for baseline benevolence shouldn’t systematically move the treatment estimate; when the covariate strongly predicts the outcome, it can substantially improve precision.

Activity 2.4 — Adjusted means: the CI and the PI

In Video 2 you plugged the average baseline into the model’s equation and got the two adjusted means by hand: 3.207 for Control and 3.710 for Treatment. predict() does that arithmetic for you and adds an interval. It works the way it did in M10: hand it the fitted model and a small tibble describing the participants you want predictions for, one row per profile. Here the two profiles are a Control participant and a Treatment participant, both at the average baseline benevolence.

Two kinds of interval answer two questions. A confidence interval (CI) describes uncertainty about the predicted mean T2 at each profile. A prediction interval (PI) describes where one new participant’s T2 might land. They come from the same call with one argument changed, exactly as in M10.

Your task: fill three blanks — the value that puts centered baseline benevolence at its average, and the kind of interval in each call.

Centered baseline benevolence is 0 for someone at the average. The interval for the mean response at a profile is a confidence interval, interval = "confidence"; the interval for one new observation is a prediction interval, interval = "prediction".

predict() returns three columns, one row per profile: fit, the prediction, and lwr and upr, the lower and upper ends of the interval. bind_cols() attaches them to profiles, so each prediction sits next to the profile it describes. (We simply wrote the labels "Control" and "Treatment": once m1 is fit, it remembers treat’s categories, so predict() only needs the label.)

Compare the two tibbles. fit is the same in both: 3.207 is the intercept \(b_0\) and 3.710 is \(b_0 + b_2\), the numbers from Video 2. Only lwr and upr differ, and by a lot. Each CI is about 0.07 of a point wide; each PI is about 3.22 points wide, roughly 45 times as wide, because it adds the person-to-person scatter (\(\sigma\) = 0.82) to the uncertainty about the mean.

Activity 2.5 — From predict() to ggplot

Turning predict() output into a graph takes one idea: every column the graph needs is already in the tibble. fit says where the dot goes; lwr and upr say where its interval starts and stops. So map y = fit, ymin = lwr, and ymax = upr.

Your task: fill three blanks — the three mappings in the CI graph. The second graph is filled in for you: it draws the PIs and the CIs together on one y-axis, using the same three columns.

In adj_means, the prediction is the column fit, and the interval’s ends are lwr and upr.

Look at the two y-axes first. The CI graph spans well under one benevolence point, so its intervals look generous; on the second graph’s shared scale, each CI shrinks to a dark sliver inside a pale PI about three points tall. The two intervals answer two questions. The CI describes uncertainty about the predicted mean T2 for each profile (each arm at the average baseline). Those intervals are short, so each profile’s mean is estimated precisely; how precisely the Treatment − Control difference is estimated comes straight from the treatment coefficient’s own CI, [0.452, 0.554]. The PI answers a question about one person: where might a new participant’s score land? The two PIs overlap almost entirely, so knowing someone’s arm tells you much less about their own score than the CIs alone suggest. That’s the lesson from M09’s boulder study: readers shown only CIs judge an effect bigger than readers shown PIs. Report the interval that matches your question. (The Treatment PI’s upper end, 5.32, passes 5, the top of the scale; M12 looks at why.)

Quick Check

Answer each question — you’ll see green (correct) or pink (incorrect) feedback as you type.

1. When you center a continuous predictor (subtract its mean from every value) before fitting a regression, what changes?

2. True or False: In a randomized trial, adjusting for a strong baseline covariate is a clever way to shift the treatment-effect estimate in the direction that supports your hypothesis.

3. Why does adjusting for baseline benevolence shrink the standard error on the treat coefficient compared to a treat-only model?

4. True or False: In the baseline-and-treatment model, the baseline-benevolence row and the treatment row of tidy() use the same critical value, \(t^*\).

5. In that model, what does the p-value in the treatTreatment row test?

6. A counselor asks: “If my next client is given the workbook, what range should I expect their T2 benevolence to fall in?” Which interval answers that question?


Video 3 — Respecting the Design: Site and the Treatment Estimate

Video looks blurry? Click the gear icon at the lower right of the player and set Quality to 1080p, or the highest option offered.

What to listen for:

  • Why add site: it was the trial’s randomization stratum, not a confounder (a variable that affects both who ends up in the treatment group and the outcome, so it can create a difference between the groups on its own) — the arms are balanced across sites by design. It belongs in the model because the model should reflect how the study was run, and because sites differ in their average scores. Each site gets its own adjusted mean level, relative to Colombia, the reference site (R builds five indicators for the other sites). With only six sites, these indicators are the simple, transparent choice; a multilevel model, which treats site as a grouping level and asks how much sites vary, is the alternative PSY 653 covers
  • Reading a site coefficient: Hong Kong’s, about −0.32, compares Hong Kong with Colombia for participants in the same arm with the same baseline
  • Why some coefficients move when a predictor enters: part of a raw gap can be carried over from where people started. The arms barely differed in baseline benevolence (randomization), so treat barely moves; the sites did differ, so the site gaps shrink toward zero
  • The nested-model F-test (anova(m1, m2)) tests whether a block of new predictors jointly adds information; ΔR² is how much it adds
  • The Table 2 fallacy: control-variable coefficients are not automatically causal effects. For site there’s a second issue: site comes before baseline, so baseline benevolence sits partly on the path from site to T2, and adjusting for it removes the part of each site gap that runs through baseline. Random assignment provides a basis for interpreting the treatment contrast causally, subject here to additional missing-follow-up assumptions
  • A category has no average, so a prediction from this model has to name a site. At every site the model-implied Treatment − Control difference is the same, the treat coefficient: site moves both predictions up or down together (the Hong Kong term cancels). The model assumes that; checking whether the effect really varies by site takes a treatment × site interaction (PSY 653)
  • The treatment estimate is an assignment-based, complete-case treatment contrast: participants are analyzed in the arm they were assigned to, but only those with a T2 score are included, so it is not the full-sample intention-to-treat analysis the published paper used. Principled missing-data methods, which can account for loss to follow-up more fully, come in more advanced courses

Activity 3.1 — Add site to the model

Add site to the baseline-and-treatment model’s predictors. R will automatically build the five site dummies (one per non-reference site). The reference site is Colombia: the page’s setup code already ran mutate(site = fct_relevel(site, "Colombia")), the same step the Module shows, so each site coefficient compares that site with Colombia.

Your task: fill one blank — the site predictor.

Site is the recruitment-site variable in REACH (six sites in five countries; Ukraine has two); the column is just called site.

You should see eight rows — the intercept, the slope on baseline benevolence, the treat slope, and five site dummies. Note that the treat coefficient is about 0.500, essentially unchanged from 0.503 without site, as we’d expect under random assignment within site. This site-indicator model follows the same logic of site adjustment the published REACH analysis used. (The published analysis differs in its details — a different primary outcome (the full unforgiveness score, not the benevolence subscale), a different estimator, no baseline covariate in the primary model, and multiple imputation rather than complete cases — so treat ours as the same idea, not the same specification.)

The glance() row summarizes the model as a whole. R² is .400: the model now accounts for about 40% of the variation in T2 benevolence. The two df columns follow the rule from Video 2. df = 7 is the number of slopes, k — 1 for baseline benevolence, 1 for treatment, and 5 for the site dummies — and df.residual = 3,957 is n − k − 1 = 3,965 − 7 − 1.

One caution before the next activity. This \(F_\text{obs}\) (about 376) tests whether all 7 slopes are 0 at once, so its enormous size mostly reflects baseline benevolence, which we already knew predicts T2. It says nothing specific about site. To ask whether the site block adds anything beyond baseline benevolence and treatment, you need the nested-model F-test in Activity 3.2.

Activity 3.2 — Nested-model F-test

Test whether the site block — the five new dummies — jointly adds information. The R function for nested-model comparison is anova().

Your task: fill two blanks — the smaller (nested) model and the larger model.

The smaller model is the one without site (m1); the larger is the one with site (m2).

\(F_\text{obs}\) is about 14.1 with 5 numerator df (the five site dummies) and 3,957 denominator df, p < .001. The data provide evidence against the hypothesis that all five site coefficients are jointly zero in the population. The F-test assesses evidence of added predictive information; ΔR² tells us how much site adds. Site belongs in the classroom model because of the design either way. The corresponding change in R² when site is added — about .011 — is the incremental R² for the site block (the block-level analog of the single-predictor \(sr^2\)).

Activity 3.3 — One site, then every site

Video 3 asked what the treatment estimate means once site is in the model. Site is a category, so it has no average, and any prediction from m2 has to name a site. Start with one site, Colombia, then do all six at once.

Your task: fill three blanks — the site for Step 1, and the two columns to subtract in Step 2.

The reference site is spelled "Colombia". After pivot_wider(), the two arms are columns named Treatment and Control; the Treatment − Control difference subtracts the second from the first.

Step 1. For a Colombia participant at the average baseline, the model predicts 3.336 under Control and 3.836 under Treatment. Colombia is the reference site, so the Control prediction is the intercept. Those numbers describe one site, about 27% of the sample.

Step 2. expand_grid() builds every combination of the values you give it: 6 sites × 2 arms = 12 profiles, all at the average baseline. pivot_wider(), from M04, turns the two arms into side-by-side columns. Read down the Control and Treatment columns: the predictions differ from site to site, from 3.02 to 3.34 under Control. Now read down difference: 0.500 in every row, the treat coefficient from Activity 3.1. As Video 3 showed with Hong Kong, the site term appears in both predictions and cancels, so only \(b_2\) is left.

So the treatment estimate from this model needs no extra calculation: it’s the treat coefficient, about 0.500 benevolence points, 95% CI [0.449, 0.550]. It is an assignment-based, complete-case contrast: everyone is analyzed in the arm they were assigned to, but only participants with a T2 score are included, so a causal reading also depends on assumptions about the missing follow-ups. The same gap at every site is something the model assumes, not something it found; whether the workbook works equally well everywhere is a moderation question for PSY 653.

The Module’s Going further box shows a second way to get adjusted means from this model: pick no site, let everyone keep their own, and average. Make two copies of analysis with treat overwritten, predict for every row, and take the mean:

These averaged means describe the sample’s mix of sites rather than Colombia, so they sit lower than Colombia’s. Their difference is the same 0.500.

Quick Check

Answer each question — you’ll see green (correct) or pink (incorrect) feedback as you type.

1. Why include site in the model, given that random assignment within site already balances treatment across sites?

2. The Table 2 fallacy is the mistake of:

3. When baseline benevolence is added to a model with treat and site, the treat coefficient barely changes, but every site coefficient shrinks toward zero. Why?

4. True or False: In the model with site, the predicted Treatment − Control difference is larger at Colombia, the highest-scoring site, than at Hong Kong.

5. In the Module’s adjusted-means-by-site figure (also shown in Video 3), the Treatment − Control gap is identical in all six panels. What does that show?


Five things to carry into lecture

  1. In the additive linear models we fitted, each slope is an adjusted comparison: the predicted difference between people one unit apart on that predictor who match on the others. Each slope automatically accounts for (“partials out”) the others; for treatment here, the outcome-only residual gap from Video 1 closely previews that adjusted comparison, and the MLR estimates it exactly.
  2. Baseline adjustment improved precision in our REACH analysis. Under randomization, baseline benevolence is balanced across arms in expectation; here, because it strongly predicts T2, adding it reduced the treatment coefficient’s standard error with little change in its estimate. This need not happen to the same degree in every study.
  3. A CI and a PI answer different questions. predict() returns both: interval = "confidence" for the predicted mean at a profile (an arm at the average baseline), and interval = "prediction" for where one new participant’s score might land. The PI is far wider because it adds person-to-person scatter, and showing only CIs can make an effect look bigger than it is for any one person.
  4. Define the estimand before choosing how to report the model. In our additive models the treatment coefficient directly estimates the adjusted difference, at any baseline and at any site. predict() gives the adjusted means for profiles you name; with a categorical predictor like site, the profile has to name a category.
  5. Control-variable coefficients are not automatically causal effects. That’s the Table 2 fallacy. Report the focal estimand (our assignment-based, complete-case treatment contrast) with its uncertainty, and recognize that a causal interpretation still depends on assumptions about missing follow-up outcomes.

Monday’s lecture builds on these steps. Bring questions — especially about the Table 2 fallacy, a common stumbling point in published work.