Multiple regression (>1 predictor)

Description and model specification

Multiple linear regression involves looking at one continuous outcome (i.e., DV), with two or more independent variables (i.e., IVs).

A multiple linear regression model takes the following form:

\[ y = \beta_0 + (\beta_1 \cdot x_{1}) + (\beta_2 \cdot x_{2}) + .... + (\beta_j \cdot x_{j}) + \epsilon \quad \text{where} \quad \epsilon \sim N(0, \sigma) \text{ independently} \]

$$
y = \beta_0 + (\beta_1 \cdot x_{1}) + (\beta_2 \cdot x_{2}) + .... + (\beta_j \cdot x_{j}) + \epsilon
\quad \text{where} \quad \epsilon \sim N(0, \sigma) \text{ independently}
$$

In R:

Multiple and simple linear regression follow the same structure within the lm() function - the logic scales up to however many predictor variables we want to include in our model. You simply add (using the + sign) more independent variables. For example, if we wanted to build a multiple linear regression that included three independent variables, we could fit one of the following via the lm() function:

model_name <- lm(DV ~ IV1 + IV2 + IV3, data = data_name)

Interpretation of coefficients

You’ll hear a lot of different ways that people explain multiple regression coefficients.

For the model \(y = \beta_0 + \beta_1 \cdot x_1 + \beta_2 \cdot x_2 + \epsilon\), the estimate \(\hat \beta_1\) will often be reported as:

“the increase in \(y\) for a one unit increase in \(x_1\) when…”

“holding the effect of \(x_2\) constant.”
“controlling for differences in \(x_2\).”
“partialling out the effects of \(x_2\).”
“holding \(x_2\) equal.”
“accounting for effects of \(x_2\).”

For models with 3+ predictors, just like building the model in R, the logic of the above simply extends.

For example “the increase in [outcome] for a one unit increase in [predictor] when…”

“holding [other predictors] constant.”
“accounting for [other predictors].”
“controlling for differences in [other predictors].”
“partialling out the effects of [other predictors].”
“holding [other predictors] equal.”
“accounting for effects of [other predictors].”

Example

Imagine that you were tasked to investigate whether recall accuracy was associated with recall confidence and age. You have been provided with data from twenty participants who studied passages of text (c500 words long), and were tested a week later. The testing phase presented participants with 100 statements about the text. They had to answer whether each statement was true or false, as well as rate their confidence in each answer (on a sliding scale from 0 to 100). The dataset contains, for each participant, the percentage of items correctly answered, their age (in years), and their average confidence rating.

The data are available at https://uoepsy.github.io/data/recalldata.csv

recalldata <- read_csv('https://uoepsy.github.io/data/recalldata.csv')

Research Question

Is recall accuracy associated with recall confidence and age?

Visualise data

p_acc <- recalldata |>
  ggplot(aes(x = recall_accuracy)) +
  geom_histogram(binwidth=5) +
  xlab('Recall\naccuracy') +
  ylab('Count')

p_con <- recalldata |>
  ggplot(aes(x = recall_confidence)) +
  geom_histogram(binwidth=5) +
  xlab('Recall\nconfidence') +
  ylab('')

p_age <- recalldata |>
  ggplot(aes(x = age)) +
  geom_histogram(binwidth=5) +
  xlab('Age') +
  ylab('')

p_acc + p_con + p_age

Figure 1: Association between Recall Accuracy, Recall Confidence, and Age

Model Specification

\[ \text{Recall Accuracy} = \beta_0 + (\beta_1 \cdot \text{Recall Confidence}) + (\beta_2 \cdot \text{Age}) + \epsilon \\ \]

$$
\text{Recall Accuracy} = \beta_0 + (\beta_1 \cdot \text{Recall Confidence}) + (\beta_2 \cdot \text{Age}) + \epsilon \\
$$

Hypothesis Specification

\(H_0:\) All \(\beta_j = 0\) (for \(j = 1, 2\))

$H_0:$ All $\beta_j = 0$ (for $j = 1, 2$)

There is no association between recall accuracy and recall confidence and/or age.

\(H_1:\) At least one \(\beta_j \neq 0\) (for \(j = 1, 2\))

$H_1:$ At least one $\beta_j \neq  0$ (for $j = 1, 2$)

There is an association between recall accuracy and recall confidence and/or age.

Fitting model

recall_multi <- lm(
  recall_accuracy ~ recall_confidence + age,
  data = recalldata
)

Interpreting results

summary(recall_multi)


Call:
lm(formula = recall_accuracy ~ recall_confidence + age, data = recalldata)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.194  -5.175  -0.553   2.593  18.681 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         36.160     12.841    2.82  0.01190 *  
recall_confidence    0.896      0.191    4.68  0.00021 ***
age                 -0.339      0.153   -2.21  0.04098 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.67 on 17 degrees of freedom
Multiple R-squared:  0.603, Adjusted R-squared:  0.557 
F-statistic: 12.9 on 2 and 17 DF,  p-value: 0.000387

\(\beta_0\) = (Intercept) = 36.16

An individual aged 0 years with no recall confidence was expected to have a recall accuracy of \(36.16\).

Note: the intercept isn’t very useful here at all. It estimates the accuracy for a newborn (who wouldn’t be able to complete the task!).

\(\beta_1\) = recall_confidence = 0.9

Holding age constant, increasing one unit in recall confidence is associated with a \(0.9\) percentage point increase in recall accuracy.
This estimate is significantly different from zero \((p < .001)\).

\(\beta_2\) = age = -0.34

Holding recall confidence constant, increasing one year of age is associated with a \(-0.34\) percentage point decrease in recall accuracy.
This estimate is significantly different from zero \((p = .041)\).

Visualising the model

We can use the function plot_model() from sjPlot.

p_con_mod <- plot_model(recall_multi,
           type = "eff",
           terms = "recall_confidence",
           show.data = TRUE)

p_age_mod <- plot_model(recall_multi,
           type = "eff",
           terms = "age",
           show.data = TRUE) +
  ylab('')

p_con_mod + p_age_mod

Figure 2: Association between Recall Accuracy, Recall Confidence, and Age