04: Model comparison

This week, you’ll practice interpreting the F-test for model fit, as well as comparing different sorts of models to one another using incremental F-tests and AIC/BIC.

NoteGet set up

1. Open RStudio.
2. Create a new R Markdown (.Rmd) file for this week’s exercises.
3. Save it somewhere you can find it again.
4. Give it a clear name (for example, dapr2_lab04.Rmd).
5. In the first code chunk, load the packages you’ll need this week (and install them if you don’t have them already):
   • tidyverse

---

Data dictionary:

variable	description
age	Age in years of respondent
outdoor_time	Self report estimated number of hours per week spent outdoors
social_int	Self report estimated number of social interactions per week (both online and in-person)
routine	Binary 1=Yes/0=No response to the question 'Do you follow a daily routine throughout the week?'
wellbeing	Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS), a self-report measure of mental health and wellbeing. The scale is scored by summing responses to each item, with items answered on a 1 to 5 Likert scale. The minimum scale score is 14 and the maximum is 70
location	Location of primary residence (City, Suburb, Rural)
steps_k	Average weekly number of steps in thousands (as given by activity tracker if available)

NoteMore detail about this dataset

From the Edinburgh & Lothians, 100 city/suburb residences and 100 rural residences were chosen at random and contacted to participate in the study. The Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS) was used to measure mental health and wellbeing.

Participants filled out a questionnaire including items concerning: estimated average number of hours spent outdoors each week, estimated average number of social interactions each week (whether on-line or in-person), whether a daily routine is followed (yes/no). For those respondents who had an activity tracker app or smart watch, they were asked to provide their average weekly number of steps.

Test a model’s fit

Question 1

Once again, read in the data from https://uoepsy.github.io/data/wellbeing_rural.csv and store it in a variable named mwdata (“mw” stands for “mental wellbeing”).

For this week’s lab, we’ll need to start by deleting all observations that contain any NA values in the following columns:

• wellbeing
• age
• outdoor_time
• social_int
• steps_k

Using the tidyverse function drop_na(), remove rows containing NAs, and then assign the resulting dataframe back to the variable name mwdata.

You’ll know you’ve done it right if mwdata now contains 134 observations, not 200.

TipCode hint

Replace the ...s with each actual variable name.

mwdata <- mwdata |>
  drop_na(..., ..., ..., ..., ...)

Solution

Solution 1.

mwdata <- read_csv("https://uoepsy.github.io/data/wellbeing_rural.csv")

mwdata <- mwdata |>
  drop_na(wellbeing, age, outdoor_time, social_int, steps_k)

nrow(mwdata)

[1] 134

Question 2

This linear model should be familiar to you by now:

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{outdoor\_time}) + (\beta_2 \cdot \text{social\_int}) + \epsilon \]

Using the function lm(), fit this linear model and name the result m0.

🗂️ See Multiple regression > Fitting model flash card.

Solution

Solution 2.

m0 <- lm(wellbeing ~ outdoor_time + social_int, data = mwdata)

Question 3

Now we want to know: is m0 significantly better at predicting wellbeing than a model with no predictors, i.e., a null model?

This question is addressed by the F-test reported at the very bottom of the model summary for m0.

Write a sentence to report the answer using APA notation.

🗂️ See Assessing model fit > F-ratio > Report in APA format flash card.

Solution

Solution 3.

summary(m0)

Call:
lm(formula = wellbeing ~ outdoor_time + social_int, data = mwdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.4504  -2.8922  -0.5264   3.2577  11.1911 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  30.11090    1.81348  16.604  < 2e-16 ***
outdoor_time  0.22915    0.06271   3.654 0.000372 ***
social_int    0.18660    0.10913   1.710 0.089646 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.203 on 131 degrees of freedom
Multiple R-squared:  0.1082,    Adjusted R-squared:  0.09461 
F-statistic: 7.949 on 2 and 131 DF,  p-value: 0.0005517

The linear model which includes the predictors outdoor time and social interaction explained significantly more variance in wellbeing score than a null model (\(F(2, 131) = 7.95\), \(p < .001\)).

Compare nested models

Question 4

New RQ: Is wellbeing significantly predicted by number of steps (steps_k) and hours outdoors (outdoor_time), over and above the impact of years of age (age) and number of social interactions (social_int)?

To address this RQ, you’ll need to compare two models.

The outcome variable for both models is wellbeing. Write down two lists (the order does not matter):

1. the predictors in one of the models
2. the predictors in the other model

Bonus challenge: Write out the formulae for these two models using mathematical model notation (that is, the format with the \(\beta\)s).

Solution

Solution 4. One model will predict wellbeing based on

• age and
• social_int.

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{age}) + (\beta_2 \cdot \text{social\_int}) + \epsilon \]

$$
\text{wellbeing} = \beta_0 + 
(\beta_1 \cdot \text{age}) +
(\beta_2 \cdot \text{social\_int}) +
\epsilon
$$

The other model will predict wellbeing based on

• age,
• social_int,
• steps_k, and
• outdoor_time.

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{age}) + (\beta_2 \cdot \text{social\_int}) + (\beta_3 \cdot \text{steps\_k}) + (\beta_4 \cdot \text{outdoor\_time}) + \epsilon \]

$$
\text{wellbeing} = \beta_0 + 
(\beta_1 \cdot \text{age}) +
(\beta_2 \cdot \text{social\_int}) +
(\beta_3 \cdot \text{steps\_k}) +
(\beta_4 \cdot \text{outdoor\_time}) +
\epsilon
$$

Question 5

Using lm(), fit the two models you identified in Q4.

Name the model with fewer predictors m1 and the model with more predictors m2.

🗂️ See Multiple regression > Fitting model flash card.

Solution

Solution 5.

m1 <- lm(wellbeing ~ age + social_int, data = mwdata)
m2 <- lm(wellbeing ~ age + social_int + steps_k + outdoor_time, data = mwdata)

Note: The predictors do not need to be in this specific order. As long as they’re in the model at all, that’s what matters.

Question 6

Compare m1 and m2 using an incremental F-test.

Based on the results of the test, write a sentence that addresses the research question and reports the incremental F-test using APA notation.

🗂️ See Comparing LMs > Incremental F-test flash card.

Solution

Solution 6.

anova(m1, m2)

Analysis of Variance Table

Model 1: wellbeing ~ age + social_int
Model 2: wellbeing ~ age + social_int + steps_k + outdoor_time
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    131 3886.8                                
2    129 3494.0  2    392.83 7.2517 0.001036 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yes, wellbeing is predicted by number of steps and hours outdoors, over and above the impact of years of age and number of social interactions (\(F(2, 129) = 7.25, p = .001\)).

Compare non-nested models

Question 7

New RQ: Do number of steps (steps_k) and hours outdoors (outdoor_time) predict wellbeing better than years of age (age) and number of social interactions (social_int)?

To address this RQ, you’ll need to compare two models.

The outcome variable for both models is wellbeing. Write down two lists (the order does not matter):

1. the predictors in one of the models
2. the predictors in the other model

Bonus challenge: Write out the formulae for these two models using mathematical model notation (that is, the format with the \(\beta\)s).

Solution

Solution 7. One model will predict wellbeing as a function of

• steps_k and
• outdoor_time.

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{steps\_k}) + (\beta_2 \cdot \text{outdoor\_time}) + \epsilon \]

$$
\text{wellbeing} = \beta_0 + 
(\beta_1 \cdot \text{steps\_k}) +
(\beta_2 \cdot \text{outdoor\_time}) +
\epsilon
$$

The other model will predict wellbeing as a function of

• age and
• social_int.

\[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{age}) + (\beta_2 \cdot \text{social\_int}) + \epsilon \]

$$
\text{wellbeing} = \beta_0 + 
(\beta_1 \cdot \text{age}) +
(\beta_2 \cdot \text{social\_int}) +
\epsilon
$$

Question 8

Using lm(), fit the two models you identified in Q7.

Name one model m3 and the other m4 (it doesn’t matter which is which).

🗂️ See Multiple regression > Fitting model flash card.

Solution

Solution 8.

m3 <- lm(wellbeing ~ steps_k + outdoor_time, data = mwdata)
m4 <- lm(wellbeing ~ age + social_int, data = mwdata)

Note: The predictors do not need to be in this specific order. As long as they’re in the model at all, that’s what matters.

Question 9

Using AIC and BIC, compare these two models.

Based on the results, write a sentence that addresses the research question, and report each model’s AIC and BIC using APA notation.

🗂️ See Comparing LMs > AIC & BIC flash card.

Solution

Solution 9. For both AIC and BIC, smaller values indicate a model with better fit.

AIC(m3, m4)

   df      AIC
m3  4 829.7469
m4  4 839.5217

• Based on AIC, the better-fitting model is m3.

BIC(m3, m4)

   df      BIC
m3  4 841.3383
m4  4 851.1131

• Based on BIC, the better-fitting model is also m3.
• And additionally, since the difference is just under 10, we have strong evidence of a difference between the two models.

Yes, a model containing number of steps and hours outdoors (AIC = 829.75, BIC = 841.34) predicts wellbeing better than a model containing years of age and number of social interactions (AIC = 839.52, BIC = 851.11).

Identify nestedness of models

Question 10

For each of the following pairs of models, state whether they are nested or not. Assume all models are fitted to the same dataset mwdata.

Pair 1:

• wellbeing ~ social_int
• wellbeing ~ social_int + outdoor_time

Pair 2:

• wellbeing ~ social_int + outdoor_time
• wellbeing ~ social_int + age

Pair 3:

• wellbeing ~ social_int
• wellbeing ~ social_int + outdoor_time + age

Pair 4:

• wellbeing ~ social_int + age
• wellbeing ~ outdoor_time + age + social_int + steps_k

Pair 5:

• wellbeing ~ social_int
• wellbeing ~ outdoor_time

🗂️ See Comparing LMs > Nested vs. non-nested models flash card.

Solution

Solution 10.

• Pair 1: Nested (all variables in one model are contained in the other model).
• Pair 2: Not nested (neither model contains all the variables that the other one contains).
• Pair 3: Nested.
• Pair 4: Nested (the order of predictors doesn’t matter, just their presence/absence).
• Pair 5: Not nested.

Optional: Knit your report to PDF

Note: During the labs, please spend your time on the exercises above, not on rendering your document. This guidance is here in case you have extra time at the end.

Rmarkdown is a useful format because it can be used as a foundation for generating different kinds of documents. For the group report, you’ll be asked to render your Rmd file to PDF. Why not practice now?

First, make sure you have tinytex installed in R. This is what lets you convert things into PDF format.

Run the following code to install tinytex (you’ll only need to run it once).

install.packages("tinytex")
tinytex::install_tinytex()

Then ensure that the “yaml” (rhymes with “mammal”) header—the bit at the top of your Rmd document—looks something like this:

---
title: "this is my report title"
author: "B1234506"
date: "07/09/2025"
output: bookdown::pdf_document2
---

Then follow the instructions from the Rmd bootcamp to knit your Rmd document to PDF.

If you encounter errors, refer to the DAPR1 formatting resources / knitting checklist.

TipWhat to do if you cannot knit to PDF

If you are having issues knitting directly to PDF, try the following:

• Knit to HTML file
  
• Open your HTML in a web-browser (e.g. Chrome, Firefox)
  
• Print to PDF (Ctrl+P, then choose to save to PDF)
  
• Open file to check formatting

TipHiding code and/or output from code chunks

Review Lesson 5 of the rmd bootcamp for a detailed description/worked examples.

Hiding R Code

To not show the code of an R code chunk, and only show the output, write:

```{r, echo=FALSE}
# code goes here
```

Hiding R Output:

To show the code of an R code chunk, but hide the output, write:

```{r, results='hide'}
# code goes here
```

Hiding R code AND output

To hide both code and output of an R code chunk, write:

```{r, include=FALSE}
# code goes here
```