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_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 cities and suburbs. 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 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), 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.

Data Dictionary

The data in wellbeing_rural.csv contain seven 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
location	Location of primary residence (City, Suburb, Rural)
steps_k	Average weekly number of steps in thousands (as given by activity tracker if available)

Preview

The first six rows of the data are:

age	outdoor_time	social_int	routine	wellbeing	location	steps_k
28	12	13	1	36	rural	21.6
56	5	15	1	41	rural	12.3
25	19	11	1	35	rural	49.8
60	25	15	0	35	rural	NA
19	9	18	1	32	rural	48.1
34	18	13	1	34	rural	67.3

Setup

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

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

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.

If we fit the interaction x1:x2, we almost always want to also fit the separate effects x1 and x2:

“Except in special circumstances, a model including a product term for interaction between two explanatory variables should also include terms with each of the explanatory variables individually, even though their coefficients may not be significantly different from zero. Following this rule avoids the logical inconsistency of saying that the effect of $X_1$ depends on the level of $X_2$ but that there is no effect of $X_1$.”
— Ramsey and Schafer (2012)

Solution

Question 2

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

Note that the “location” variable currently has three levels (Rural/Suburb/City). In order to address the research question, we only want two (Rural/Not Rural) locations - you will need to fix this.

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

Hint

One way to do this would be to use ifelse() to define a variable which takes one value (“Rural”) if the observation meets from some condition, or another value (“Not Rural”) if it does not.

Type ?ifelse in the console if you want to see the help function. You can use it to add a new variable either inside mutate(), or using data$new_variable_name <- ifelse(test, x, y) syntax.

Solution

Descriptive Statistics & Visualisations

Question 3

Provide a table of descriptive statistics and visualise your data.

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

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.
To include facets, within your ggplot() argument you will need to specify + facet_wrap() in order to produce facets for each location. It would also be useful to specify geom_smooth(method="lm").
To include different colours for each location, you will need to specify colour = isRural within your ggplot() argument, as well as + scale_colour_discrete() to state the name = of your legend and the labels = associated with each coloured line.

Solution

Numeric
Visual

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

mwdata %>%
    select(isRural, wellbeing, social_int) %>%
    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(mwdata, mwdata$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.0  70.0  52.0  0.13
outdoor_time    2 100 18.72  6.91   18.0   18.61  7.41  6.0  34.0  28.0  0.08
social_int      3 100 11.67  3.95   12.0   11.51  3.71  3.0  24.0  21.0  0.37
routine         4 100  0.54  0.50    1.0    0.55  0.00  0.0   1.0   1.0 -0.16
wellbeing       5 100 38.57  5.66   39.0   38.36  5.93 26.0  59.0  33.0  0.44
location*       6 100  1.43  0.50    1.0    1.41  0.00  1.0   2.0   1.0  0.28
steps_k         7  68 51.03 26.14   49.6   50.56 31.58  2.9 111.3 108.4  0.15
isRural*        8 100  1.00  0.00    1.0    1.00  0.00  1.0   1.0   0.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
location*       -1.94 0.05
steps_k         -0.94 3.17
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.00   41.38 17.79  18  69.0  51.0  0.09
outdoor_time    2 100 17.79  7.29  18.00   17.68  7.41   1  35.0  34.0  0.06
social_int      3 100 12.46  4.08  12.00   12.41  4.45   4  21.0  17.0  0.05
routine         4 100  0.59  0.49   1.00    0.61  0.00   0   1.0   1.0 -0.36
wellbeing       5 100 34.02  3.99  34.00   34.08  2.97  22  45.0  23.0 -0.07
location*       6 100  1.00  0.00   1.00    1.00  0.00   1   1.0   0.0   NaN
steps_k         7  66 38.63 26.64  35.15   37.00 30.62   0 110.6 110.6  0.45
isRural*        8 100  2.00  0.00   2.00    2.00  0.00   2   2.0   0.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
location*         NaN 0.00
steps_k         -0.59 3.28
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):

mwdata %>% 
  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 = mwdata, 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 = mwdata, 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 = mwdata, 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 your model using lm(), and assign it as an object with the name “rural_mod”.

Hint

When fitting a regression model in R with two explanatory variables A and B, and their interaction, these three are equivalent:

y ~ A + B + A:B
y ~ A + B + A*B
y ~ A*B

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.

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. We could either create a new variable, named say mc_social_int which is social_int - mean(social_int) (which would therefore make everyone who had 12.06 social interactions per week be given a value of 0 in the variable), or we could do it in the model:

Creating Mean Centred Variable
Mean Centre in Model

mwdata <-
 mwdata %>%
  mutate(
   mc_social_int = social_int - mean(social_int)
    )

rural_mod_mc <- lm(wellbeing ~ I(social_int-12.06) * isRural, data = mwdata)

What does the I do?

When we specify the I() notation in the formula syntax of lm() in R, the notation simply means ‘as is’. In other words, in our example, we are specifying that we want to include as a predictor variable the product of social_int - 12.06 in our model.

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 = mwdata)

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 = mwdata)

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 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 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

In terms of specification, it might be useful to look up the helper function via ?probe_interaction. As a quick guide:

model = name model to be used
pred = The predictor variable
modx = The moderator variable
interval = If we say TRUE, then confidence/prediction intervals will be plotted around the line

Make sure to give your plot informative titles/labels. You, for example, likely want to give your plot:

a clear and concise title (specify main.title =)
axis labels with units or scale included (specify x.label = and y.label =)
a legend title (specify legend.main =)

Solution

Writing Up & Presenting Results

Question 9

Provide key model results (from “rural_mod1”) in a formatted table.

Hint

Use tab_model() from the sjPlot package.

You can rename your DV and IV labels by specifying dv.labels and pred.labels. To do so, specify your variable name on the left, and what you would like this to be named in the table on the right.

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 (from “rural_mod1”) in the context of the research question and report your model in full.

Make reference to the interaction plot and regression table.

“The best method of communicating findings about the presence of a significant interaction may be to present a table or graph of the estimated means at various combinations of the interacting variables.”
— Ramsey and Schafer (2012)

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

References

Ramsey, Fred, and Daniel Schafer. 2012. The Statistical Sleuth: A Course in Methods of Data Analysis. Cengage Learning.