Interactions I: Num x Cat

Learning Objectives

At the end of this lab, you will:

Understand the concept of an interaction
Be able to interpret the meaning of a numeric \(\times\) categorical interaction
Be able to visualize and probe interactions

What You Need

Be up to date with lectures
Have completed all labs from Semester 1

Required R Packages

Remember to load all packages within a code chunk at the start of your RMarkdown file using library(). If you do not have a package and need to install, do so within the console using install.packages(" "). For further guidance on installing/updating packages, see Section C here.

For this lab, you will need to load the following package(s):

tidyverse
kableExtra
psych
sjPlot
patchwork
sandwich
interactions

Presenting Results

All results should be presented following APA guidelines.If you need a reminder on how to hide code, format tables/plots, etc., make sure to review the rmd bootcamp.

The example write-up sections included as part of the solutions are not perfect - they instead should give you a good example of what information you should include and how to structure this. Note that you must not copy any of the write-ups included below for future reports - if you do, you will be committing plagiarism, and this type of academic misconduct is taken very seriously by the University. You can find out more here.

Lab Data

You can download the data required for this lab here or read it in via this link https://uoepsy.github.io/data/wellbeing_location_rural.csv

Study Overview

Research Question

Does the association between wellbeing and the number of social interactions differ between rural and non-rural residents?

Wellbeing/Rurality data codebook.

Background

Researchers were specifically interested in how the number of social interactions influenced mental health and wellbeing differently for those living in rural communities compared to those in non rural locations. They wanted to assess whether the effect of social interactions on wellbeing was moderated by (depends upon) whether or not a person lived in a rural area.

Description

From the Edinburgh & Lothians, 100 non-rural (e.g., 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 well-being.

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), and whether a daily routine is followed (yes/no).

Data Dictionary

The data in wellbeing_location_rural.csv contain six attributes collected from a random sample of \(n=200\) hypothetical residents over Edinburgh & Lothians, and include:

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 well-being. 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
isRural	Location of primary residence (Rural, Not Rural)

Preview

The first six rows of the data are:

age	outdoor_time	social_int	routine	wellbeing	isRural
28	12	13	1	36	rural
56	5	15	1	41	rural
25	19	11	1	35	rural
60	25	15	0	35	rural
19	9	18	1	32	rural
34	18	13	1	34	rural

Setup

Create a new RMarkdown file
Load the required package(s)
Read the wellbeing_location_rural dataset into R, assigning it to an object named ruraldata

Solution

Exercises

Study & Analysis Plan Overview

Question 1

Provide a brief overview of the study design and data, before detailing your analysis plan to address the research question.

Hint

Give the reader some background on the context of the study
State what type of analysis you will conduct in order to address the research question
Specify the model to be fitted to address the research question (note that you will need to specify the reference level of your categorical variable)
Specify your chosen significance (\(\alpha\)) level
State your hypotheses

Much of the information required can be found in the Study Overview codebook.

The statistical models flashcards may also be useful to refer to. Specifically the interaction models flashcards and numeric x categorical example flashcards might be of most use.

Solution

Question 2

Check coding of variables (e.g., that categorical variables are coded as factors).

As specified in Q1, we want ‘not rural’ as the reference group, so make sure to specify this.

Hint

Review the specifying reference levels flashcard.

Solution

Descriptive Statistics & Visualisations

Question 3

Provide a table of descriptive statistics and visualise your data.

Remember to interpret your findings in the context of the study (i.e., comment on any observed differences among groups).

In particular:

Explore the associations among the variables included in your analysis
Produce a visualisation of the association between weekly number of social interactions and well-being, with separate facets for rural vs non-rural respondents OR with different colours for each level of the isRural variable.

Hint

Review the many ways to numerically and visually explore your data by reading over the data exploration flashcards.

For examples, see flashcards on descriptives statistics tables - categorical and numeric values examples and numeric x categorical example - visualise data, paying particular attention to the type of data that you’re working with.

More specifically:
1. For your table of descriptive statistics, both the group_by() and summarise() functions will come in handy here.

If you use the select() function and get an error along the lines of Error in select...unused arguments..., you will need to specify dplyr::select() (this just tells R which package to use the select function from).
The pairs.panels() function from the psych package can plot all variables in a dataset against one another. This will save you the time you would have spent creating individual plots, but is only useful for continuous variables.

Solution

Numeric
Visual

We can present our summary statistics for wellbeing and social interactions grouped by location in a well formatted table using kable():

ruraldata %>%
    group_by(isRural) %>%
    summarise(Wellbeing_M = mean(wellbeing),
              Wellbeing_SD = sd(wellbeing),
              SocialInt_M = mean(social_int),
              SocialInt_SD = sd(social_int)) %>%
    kable(caption = "Wellbeing, Social Interactions, and Location Descriptive Statistics", align = "c", digits = 2, booktabs = TRUE) %>%
    kable_styling(full_width = FALSE)

Table 1: Wellbeing, Social Interactions, and Location Descriptive Statistics
isRural	Wellbeing_M	Wellbeing_SD	SocialInt_M	SocialInt_SD
not rural	38.57	5.66	11.67	3.95
rural	34.02	3.99	12.46	4.08

A quick non-APA style table - i.e., not for use in reports

For a quick overview of several summary statistics at once for all columns of our dataframe, but separately for each level of location (i.e.,rural and non-rural), you can use describeBy() from the psych package:

describeBy(ruraldata, ruraldata$isRural)


 Descriptive statistics by group 
group: not rural
             vars   n  mean    sd median trimmed   mad min max range  skew
age             1 100 43.01 14.95   41.5   42.74 18.53  18  70    52  0.13
outdoor_time    2 100 18.72  6.91   18.0   18.61  7.41   6  34    28  0.08
social_int      3 100 11.67  3.95   12.0   11.51  3.71   3  24    21  0.37
routine         4 100  0.54  0.50    1.0    0.55  0.00   0   1     1 -0.16
wellbeing       5 100 38.57  5.66   39.0   38.36  5.93  26  59    33  0.44
isRural*        6 100  1.00  0.00    1.0    1.00  0.00   1   1     0   NaN
             kurtosis   se
age             -1.13 1.49
outdoor_time    -0.87 0.69
social_int       0.13 0.39
routine         -1.99 0.05
wellbeing        0.38 0.57
isRural*          NaN 0.00
------------------------------------------------------------ 
group: rural
             vars   n  mean    sd median trimmed   mad min max range  skew
age             1 100 41.59 14.85     42   41.38 17.79  18  69    51  0.09
outdoor_time    2 100 17.79  7.29     18   17.68  7.41   1  35    34  0.06
social_int      3 100 12.46  4.08     12   12.41  4.45   4  21    17  0.05
routine         4 100  0.59  0.49      1    0.61  0.00   0   1     1 -0.36
wellbeing       5 100 34.02  3.99     34   34.08  2.97  22  45    23 -0.07
isRural*        6 100  2.00  0.00      2    2.00  0.00   2   2     0   NaN
             kurtosis   se
age             -1.18 1.49
outdoor_time    -0.49 0.73
social_int      -0.82 0.41
routine         -1.89 0.05
wellbeing        0.40 0.40
isRural*          NaN 0.00

Let’s first plot the continuous variables included within our model (note that we could use this for the whole dataset, but we don’t want to include irrelevant / non-continuous variables):

ruraldata %>% 
  select(wellbeing, social_int) %>%
  pairs.panels()

Wellbeing and social interactions appear to follow unimodal distributions. There was a weak, positive association between wellbeing and social interactions \((r = .24)\).

Now lets look at wellbeing scores by location:

ggplot(data = ruraldata, aes(x = isRural, y = wellbeing)) +
  geom_boxplot() + 
  labs(x = "Location", y = "Wellbeing (WEMWBS Scores)")

Those in rural locations appear to have lower wellbeing scores in comparison to those in non-rural locations.

Next, lets produce our plots with a facet for rural vs non-rural residents:

ggplot(data = ruraldata, aes(x = social_int, y = wellbeing)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  facet_wrap(~isRural, labeller = "label_both") + 
  labs(x = "Social Interactions (number per week)", y = "Wellbeing (WEMWBS Scores)")

Or instead of facets, we could use different colors for each location (rural vs non-rural):

ggplot(data = ruraldata, aes(x = social_int, y = wellbeing, colour = isRural)) +
  geom_point() + 
  geom_smooth(method="lm", se=FALSE) +
    scale_colour_discrete(
    name ="Location",
    labels=c("Not Rural", "Rural")) + 
    labs(x = "Social Interactions (number per week)", y = "Wellbeing (WEMWBS Scores)")

Those in non-rural locations appear to have higher wellbeing scores across almost all levels of social interactions. The slopes appear to be different for each location, where the greatest difference in wellbeing scores by location is most visible the highest number of social interactions. This suggests that there may be an interaction.

How do we know there might be an interaction?

The lines in the two plots above do not run in parallel - this suggested the presence of an interaction. Specifically in our example, the non-parallel lines suggested an interaction effect based on location, as the number of social interactions did not appear to have the same influence on rural and non-rural residents’ wellbeing scores.

However, the only way we can determine whether there is actually an interaction is by including an interaction term in our model, and testing this.

Model Fitting & Interpretation

Question 4

Fit the specified model, and assign it the name “rural_mod”.

Hint

We can fit interaction models using the lm() function.

For an overview, see the interaction models flashcards.

For an example, review the interaction models > numeric x categorical example > model building flashcards.

Solution

Question 5

Look at the parameter estimates from your model, and write a description of what each one corresponds to on the plot shown in Figure 1 (it may help to sketch out the plot yourself and annotate it, and refer to the drop down options below).

Figure 1: Multiple regression model: Wellbeing ~ Social Interactions * is Rural

Options for Mapping Parameter Estimates to Plot

Hint

Recall that we can obtain our parameter estimates using various functions such as summary(),coef(), coefficients(), etc.

For an overview of how to interpret coefficients, review the interaction models > interpreting coefficients flashcard.

For a specific example of coefficient interpretation, review the interaction models > numeric x categorical example > results interpretation flashcards.

Solution

Obtain parameter estimates:

coefficients(rural_mod)

            (Intercept)              social_int            isRuralrural 
             30.9985688               0.6487945               1.3865688 
social_int:isRuralrural 
             -0.5175856

\(\beta_0\) = (Intercept) = 31

On plot: The point at which the red line cuts the y-axis.
Interpretation: The intercept, or predicted wellbeing score when the number of social interactions per week was 0, and when location was not rural.
- A non-rural resident who had zero social interactions per week was expected to have a wellbeing score of 31.

\(\beta_1\) = social_int = 0.65

On plot: The slope (vertical increase on the y-axis associated with a 1 unit increase on the x-axis) of the red line.
Interpretation: The simple slope of social interactions (number per week) for location reference group (not rural).
- For someone who lives in a non-rural location, every 1 additional social interaction per week was associated with a 0.65 point change in their wellbeing score.

\(\beta_2\) = isRuralrural = 1.39

On plot: The vertical distance from the red to the blue line at the y-axis (where social_int = 0).
Interpretation: The simple effect of location (or the difference in wellbeing scores between rural and non rural residents) when number of social interactions was 0.
- For residents who had zero social interactions per week, living in a rural location was associated with wellbeing scores 1.39 points higher than living in non-rural locations (note that this difference was not significantly different from zero).

\(\beta_3\) = social_int:isRuralrural = -0.52

On plot: How the slope of the line differs when you move from the red to the blue line.
Interpretation: The interaction between social interactions (number per week) and location (rural/not rural). This is the estimated difference in simple slopes of social interactions for rural vs non-rural residents.
- Compared to living in non-rural locations, living in a rural location was associated with a 0.52 lesser increase in wellbeing scores for every every 1 additional social interaction per week.

Question 6

No participants in our dataset had zero hours of social interactions per week (the lowest was 3), and we’re likely not interested in differences between rural and non-rural residents who have never interacted with others.

Mean center the continuous IV(s), and re-run your model with mean centered variable(s).

Hint

There are a couple of different ways that we can re-centre. See the data transformations > centering flashcards for a recap. Note, it would be best to create a new mean-centered variable to then use within the model in this instance.

Solution

Question 7

Note any differences between the summary() output between the original (“rural_mod”) and mean centred (“rural_mod1”) models. Pay particular attention to your coefficients and their significance values. How have your coefficients changed? Why do you think these differences have been observed?

Hint

These plots illustrate the difference between the “rural_mod” and “rural_mod1” models.

Figure 2: Difference when social interactions is not vs is mean centered.
Note that the lines without SE intervals on the left plot represent predicted values below the minimum observed number of social interactions, to ensure that zero on the x-axis is visible

Solution

By comparing the summary() outputs, you should see that the coefficients for (Intercept) and isRuralrural differed between the two models, whilst the coefficients for social_int and mc_social_int were the same, as were the interaction estimates social_int:isRuralrural.

Compare Model Estimates

Original Model - rural_mod:

Call:
lm(formula = wellbeing ~ social_int * isRural, data = ruraldata)

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              30.9986     1.4284  21.702  < 2e-16 ***
social_int                0.6488     0.1160   5.593 7.42e-08 ***
isRuralrural              1.3866     2.0510   0.676  0.49981    
social_int:isRuralrural  -0.5176     0.1615  -3.206  0.00157 **

Model with mean centered social interactions - rural_mod1:

Call:
lm(formula = wellbeing ~ mc_social_int * isRural, data = ruraldata)

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 38.8263     0.4581  84.754  < 2e-16 ***
mc_social_int                0.6488     0.1160   5.593 7.42e-08 ***
isRuralrural                -4.8581     0.6478  -7.500 2.17e-12 ***
mc_social_int:isRuralrural  -0.5176     0.1615  -3.206  0.00157 **

Recall that when there is an interaction A\(\times\)B, the coefficients A and B are no longer main effects. Instead, they are conditional effects upon the other being zero.

If we have the interaction y ~ x1 + x2 + x1:x2 (where \(y\) = wellbeing; \(x_1\) = social interactions; and \(x_2\) = whether or not the respondent lives in a rural location.), then:

In our “rural_mod”, the coefficient for x2 represents the association between \(y\) and \(x_2\) for someone with a score of 0 on \(x_1\)
In our “rural_mod1”, where x1 is mean centered, this will now make the coefficient for x2 represent the association between \(y\) and \(x_2\) for someone at the average of \(x_1\) (i.e., in our current example, 12.06)

Whilst the difference in rural vs non-rural may not have been significantly different when the number of weekly social interactions is zero, there did appear to be a significant difference at the average number of social interactions (as you can see from the plot below - note that this is the same plot as in the hint). Note that we can see that the model doesn’t change, it is just extracting different information (the distance to move from the blue dot to the red dot is different):

Figure 3: Difference when social interactions is not vs is mean centered

Visualise Interaction Model

Question 8

Using the probe_interaction() function from the interactions package, visualise the interaction effects from your model.

Try to summarise the interaction effects in a short and concise sentence.

Hint

For an overview and example, review the interaction models > numeric x categorical example > model visualisation flashcards.

Solution

Writing Up & Presenting Results

Question 9

Provide key model results in a formatted table.

Hint

Use tab_model() from the sjPlot package. For a quick guide, review the tables flashcard.

Solution

Table 2: Regression Table for Wellbeing Model
	Wellbeing (WEMWBS Scores)
Predictors	Estimates	CI	p
(Intercept)	38.83	37.92 – 39.73	<0.001
Social Interactions (number per week; mean centred)	0.65	0.42 – 0.88	<0.001
Location - Rural	-4.86	-6.14 – -3.58	<0.001
Social Interactions (number per week; mean centred) * Location - Rural	-0.52	-0.84 – -0.20	0.002
Observations	200
R² / R² adjusted	0.296 / 0.285

Question 10

Interpret your results in the context of the research question and report your model in full.

Make reference to the interaction plot and regression table.

Hint

For an example of coefficient interpretation, review the interaction models > numeric x categorical example > results interpretation flashcards.

Solution

Compile Report

Knit your report to PDF, and check over your work. To do so, you should make sure:

Only the output you want your reader to see is visible (e.g., do you want to hide your code?)
Check that the tinytex package is installed
Ensure that the ‘yaml’ (bit at the very top of your document) looks something like this:

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

What 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

Hiding Code and/or Output

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

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

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

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

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

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

Tinytex

You must make sure you have tinytex installed in R so that you can “Knit” your Rmd document to a PDF file:

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

Solution