Correlation & Linear Regression: Simple and Multiple

1Objectives¶

By the end of this lesson you will be able to:

1. Compute a correlation coefficient with cor(y ~ x, data=).

2. Fit a simple linear regression model with lm(y ~ x, data=).

3. Read the coefficient table from msummary(model) and interpret the slope, the intercept, and $R^2$ in the context of the data.

4. Predict a new value with makeFun() and with predict().

5. Build and read a residual plot to check whether a straight line was a reasonable model in the first place.

6. Fit a multiple regression model with lm(y ~ x1 + x2, data=), interpret each slope “holding the other predictor(s) constant,” and compare $R^2$ to adjusted $R^2$ when a second predictor is added.

7. Recognize how a categorical predictor (a factor) enters a regression model as a set of reference-level dummy coefficients, and predict from a model that includes one.

2From “is there a relationship” to “what is the relationship”¶

L07 plotted exam_score against study_min and computed cor(exam_score ~ study_min, data = survey), getting about 0.56 — a real, positive, moderate relationship, more study time going with a higher score. Correlation answers whether (and how strongly) two numerical variables move together. Regression goes one step further and answers how: it fits a straight line that turns any study time into a predicted exam score. This lesson picks up exactly where L07 left off, same variables, same dataset — first with one predictor (simple regression, Sections 1–5), then with more than one at once (multiple regression, Section 6):

suppressMessages({library(mosaic); library(BSDA)})
set.seed(2200)
survey <- read.csv("data/survey_sim.csv")

31. Correlation, again — and what it does not tell you¶

cor(exam_score ~ study_min, data = survey)

[1] 0.5594575

cor(y ~ x, data=) is mosaic’s formula version of the correlation coefficient $r$ — one number, always between -1 and 1, summarizing the strength and direction of a linear relationship. $r = 0.56$: as weekly study time goes up, exam score tends to go up too, moderately consistently. What $r$ alone can’t do is turn “more study time” into a specific predicted score — for that you need an equation, which is exactly what lm() builds.

42. Fitting the line: lm()¶

That trend line is a least-squares regression line: the one straight line that makes the total squared vertical distance from every point to the line as small as possible. lm() (“linear model”) fits it and hands back an object holding everything about the fit:

model <- lm(exam_score ~ study_min, data = survey)
model

Call:
lm(formula = exam_score ~ study_min, data = survey)

Coefficients:
(Intercept)    study_min  
   47.34733      0.05765  

Same y ~ x, data= formula shape as every other function in this book — “model exam score as a function of study time, using survey.” Printed on its own, model gives just the two numbers that define the line:

• $b_0 = 47.347$, the intercept: the model’s predicted exam score for a student with study_min = 0. It’s a mathematical anchor point for the line, not necessarily a realistic prediction — no student in this dataset actually studied zero minutes.

• $b_1 = 0.0577$, the slope: for every additional minute of weekly study time, predicted exam score rises by about 0.0577 points, holding nothing else constant (there’s nothing else in this model to hold constant — that matters more once a model has several predictors). Scaled up, that’s about 5.8 points per extra 100 minutes of weekly study time.

53. The full picture: msummary()¶

model alone shows the line; msummary() — mosaic’s cleaner version of base R’s summary() for a fitted model — adds the coefficient table, the significance tests on each coefficient, and $R^2$:

msummary(model)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 47.347331   2.075349  22.814  < 2e-16 ***
study_min    0.057650   0.007863   7.332 3.09e-11 ***

Residual standard error: 8.157 on 118 degrees of freedom
Multiple R-squared:  0.313,	Adjusted R-squared:  0.3072 
F-statistic: 53.76 on 1 and 118 DF,  p-value: 3.093e-11

Each row of the coefficient table is its own hypothesis test — same logic as every test in L10/L11 — testing $H_0: \beta = 0$ (“this term adds nothing to the model”) against $H_a: \beta \ne 0$. The study_min row’s Pr(>|t|) of 3.09e-11 is far below $\alpha = 0.05$: reject $H_0$ — study time is a statistically significant predictor of exam score in this sample. Std. Error is the standard error of that coefficient estimate (how much the slope would jitter across repeated samples); t value is Estimate / Std. Error.

$R^2$, the Multiple R-squared line (0.313), is the proportion of the variation in exam score that this model explains using study time alone — about 31.3%. You can also pull it directly:

rsquared(model)

[1] 0.3129927

$R^2 = 0.313$ means study time accounts for a bit under a third of why exam scores vary from student to student in this sample — a real, meaningful relationship (matching the moderate $r = 0.56$ from the correlation above, since $r^2 = 0.56^2 \approx 0.31$), but far from the whole story. The other 69% comes from everything else study_min alone can’t capture: how well-rested a student was, prior background, how the questions happened to fall for them, and plain chance.

64. Predicting a new value: makeFun() and predict()¶

Once a model is fit, turning a new study time into a predicted score takes one line — two equivalent ways to do it. makeFun() (a mosaic function) turns the whole model into an ordinary R function you can call by name:

predict_exam <- makeFun(model)
predict_exam(study_min = 300)

       1 
64.64231 

A student who studies 300 minutes a week is predicted to score about 64.6. Base R’s predict() does the same job, taking a small data frame of new x values instead of named arguments — handy for predicting several values at once:

predict(model, newdata = data.frame(study_min = c(100, 300, 500)))

       1        2        3 
53.11232 64.64231 76.17230 

Both agree exactly at study_min = 300 (64.64) — makeFun() is simply a friendlier wrapper around the same predict() machinery underneath.

75. Checking the line: a residual plot¶

A residual is the gap between what actually happened and what the model predicted: $\text{residual} = \text{actual} - \text{predicted}$. resid() and fitted() pull both off a fitted model:

favstats(~resid(model))

       min        Q1    median       Q3      max         mean       sd   n
 -18.65345 -6.406816 0.5248858 5.656097 18.19894 4.706269e-16 8.122442 120
 missing
       0

The mean residual is (up to rounding) exactly 0 — least squares guarantees that, so it’s not informative on its own. What is informative is the pattern in the residuals, best seen as a plot of residuals against the model’s fitted (predicted) values:

gf_point(resid(model) ~ fitted(model), color = "#0072B2", alpha = 0.75) %>%
  gf_hline(yintercept = 0, color = "#D55E00", linetype = "dashed") %>%
  gf_labs(title = "Residual plot: exam_score ~ study_min",
          x = "Fitted (predicted) exam score", y = "Residual (actual - predicted)")

This is what a healthy residual plot looks like: points scattered randomly above and below the dashed zero line, no curve, and no funnel shape (spread staying roughly constant as the fitted values increase). That supports two of simple linear regression’s key conditions — linearity (a straight line was a reasonable choice, no missed curve) and constant variance (the model isn’t systematically worse for high or low predictions). If you instead saw a clear U-shape, a straight line would be the wrong model; if you saw the spread fanning out at one end, the constant-variance condition would be in question. Reading a residual plot is exactly this: look for any leftover pattern, because a pattern means the model missed something the straight line should have captured.

86. Multiple regression: more than one predictor¶

Section 2’s model explains exam score with study_min alone. Real questions rarely have just one plausible explanation. Multiple regression uses the same least-squares idea — minimize the total squared gap between actual and predicted $y$ — but fits a flat surface through several predictors at once instead of a single line, and estimates each predictor’s own slope after accounting for the others already in the model:

Same lm() function, same y ~ x, data= formula shape from Section 2 — just add more predictors on the right with +.

8.16.1 Two numeric predictors: does sleep time add anything?¶

sleep_hours is another numeric column in survey. Does knowing how much a student sleeps improve the prediction, on top of already knowing how much they study?

model2 <- lm(exam_score ~ study_min + sleep_hours, data = survey)
model2

Call:
lm(formula = exam_score ~ study_min + sleep_hours, data = survey)

Coefficients:
(Intercept)    study_min  sleep_hours  
   43.34753      0.05748      0.59550  

msummary(model2)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 43.347527   5.387834   8.045 7.89e-13 ***
study_min    0.057479   0.007877   7.297 3.83e-11 ***
sleep_hours  0.595501   0.740065   0.805    0.423    

Residual standard error: 8.169 on 117 degrees of freedom
Multiple R-squared:  0.3168,	Adjusted R-squared:  0.3051 
F-statistic: 27.12 on 2 and 117 DF,  p-value: 2.099e-10

With more than one predictor, read each slope as “holding the other predictor(s) constant.”

• $b_1 = 0.0575$ for study_min: holding sleep_hours fixed, each extra weekly minute of study is associated with about 0.0575 more predicted points — almost identical to Section 2’s simple-model slope (0.0577), because study_min and sleep_hours barely relate to each other in this sample, so adding sleep_hours hardly changes what study_min alone was already saying.

• $b_2 = 0.5955$ for sleep_hours: holding study_min fixed, each extra hour of sleep is associated with about 0.60 more predicted points — but Pr(>|t|) = 0.423, far above $\alpha = 0.05$. Fail to reject $H_0: \beta_{sleep\_hours} = 0$: this sample gives no evidence that sleep hours add real predictive value once study time is already in the model.

Compare $R^2$ to Section 2’s simple model:

rsquared(model)    # Section 2's model: exam_score ~ study_min
rsquared(model2)   # this model: exam_score ~ study_min + sleep_hours

[1] 0.3129927
[1] 0.3167737

Multiple R-squared barely moved, 0.313 to 0.317 — and Adjusted R-squared actually dropped, 0.3072 to 0.3051. That is not a mistake: plain $R^2$ can only go up (or stay flat) every time you add any predictor, useful or not, because the model always has one more knob to fit the training data with. Adjusted $R^2$ corrects for that by penalizing predictors that don’t pull their weight, which is exactly why it dropped here — sleep_hours isn’t earning its place in the model. Adjusted $R^2$, not plain $R^2$, is the fair number for comparing models with different counts of predictors.

Predicting works exactly as in Section 4, now supplying a value for every predictor in the model:

predict_exam2 <- makeFun(model2)
predict_exam2(study_min = 300, sleep_hours = 7)

       1 
64.75968 

predict(model2, newdata = data.frame(study_min = c(100, 300, 500),
                                      sleep_hours = c(6, 7, 8)))

       1        2        3 
52.66842 64.75968 76.85094 

8.26.2 A categorical predictor: how a factor enters lm()¶

Numeric predictors aren’t the only option — major_area is categorical, and lm() handles that automatically once the column is a factor (L11 covers factor() and why to use it):

survey$major_area <- factor(survey$major_area)
model3 <- lm(exam_score ~ study_min + major_area, data = survey)
model3

Call:
lm(formula = exam_score ~ study_min + major_area, data = survey)

Coefficients:
          (Intercept)              study_min  major_areaKinesiology  
             45.96885                0.06087                3.80364  
    major_areaNursing        major_areaOther         major_areaSTEM  
              3.57923               -3.42240               -0.57737  

A factor never gets one single coefficient — R picks its first level alphabetically (here, “Business,” since factor(major_area) wasn’t given an explicit levels = order) as the reference level, folds it into the intercept, and gives every other level its own coefficient measuring the difference from that reference, holding study_min constant. Five majors means four dummy coefficients ($k - 1$ dummies for $k$ levels) — the same “$k - 1$” you saw for major_area’s degrees of freedom in L13’s ANOVA.

msummary(model3)

                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           45.968849   2.345094  19.602  < 2e-16 ***
study_min              0.060872   0.007643   7.964 1.37e-12 ***
major_areaKinesiology  3.803636   2.137903   1.779   0.0779 .  
major_areaNursing      3.579228   2.294727   1.560   0.1216    
major_areaOther       -3.422398   2.198602  -1.557   0.1223    
major_areaSTEM        -0.577368   2.137936  -0.270   0.7876    

Residual standard error: 7.862 on 114 degrees of freedom
Multiple R-squared:  0.3834,	Adjusted R-squared:  0.3563 
F-statistic: 14.17 on 5 and 114 DF,  p-value: 8.901e-11

Read major_areaKinesiology’s row the same “holding constant” way: holding weekly study time fixed, Kinesiology majors’ predicted exam score is about 3.80 points higher than Business majors (the reference level) — though Pr(>|t|) = 0.078 is still above 0.05, so even this largest gap isn’t statistically significant at the usual threshold. major_area moves both $R^2$ (0.383, up from 0.313) and adjusted $R^2$ (0.356, up from 0.307) by noticeably more than sleep_hours did — a bigger, more genuine jump, even though no single major’s dummy coefficient clears $\alpha = 0.05$ on its own.

Predicting from a model with a factor predictor needs a value for every predictor, including the categorical one, spelled exactly as it appears in the data:

predict_exam3 <- makeFun(model3)
predict_exam3(study_min = 300, major_area = "STEM")
predict_exam3(study_min = 300, major_area = "Business")

       1 
63.65323 
      1 
64.2306 

predict(model3, newdata = data.frame(study_min = 300,
                                      major_area = c("Business", "STEM", "Nursing")))

       1        2        3 
64.23060 63.65323 67.80982 

9Summary¶

• cor(y ~ x, data=) measures the strength and direction of a linear relationship; lm(y ~ x, data=) fits the actual line, and model printed alone shows its intercept and slope.

• msummary(model) is the full picture: a coefficient table (each row its own $H_0: \beta = 0$ test), the residual standard error, and $R^2$ — rsquared(model) pulls that last number directly.

• Interpret the slope as “for each one-unit increase in $x$, predicted $y$ changes by this much”; interpret $R^2$ as “this fraction of the variation in $y$ is explained by $x$.”

• makeFun(model) turns a fitted model into a callable R function; predict(model, newdata = data.frame(x = ...)) does the same job for one or many new $x$ values — never trust either one far outside the range of $x$ the model was fit on.

• A residual plot (resid(model) vs. fitted(model)) is how you check whether a straight line was the right model at all: look for random scatter around zero, not a curve or a funnel.

• lm(y ~ x1 + x2, data=) fits a multiple regression; read each slope as “holding the other predictor(s) constant,” and compare adjusted $R^2$ (not plain $R^2$) when deciding whether an added predictor is earning its place — plain $R^2$ can only rise as you add predictors, even useless ones, but adjusted $R^2$ can fall.

• A factor predictor enters lm() as $k - 1$ dummy coefficients for its $k$ levels; R’s first level alphabetically (unless you set levels = yourself) becomes the reference level folded into the intercept, and every other level’s coefficient is its difference from that reference, holding the other predictors constant.

• The full script that generated the figures and Sections 1–5’s numbers in this lesson is committed at data/make_L12_figures.R — run it yourself to reproduce those exactly. The recap scatterplot above is the same file data/make_L07_figures.R produced in L07. Section 6’s multiple-regression numbers reuse that same seed and dataset — see the note at the end of Section 6.2 for how to reproduce them.