W3 Exercises: Nested and Crossed Structures

Psychoeducation Treatment Effects

Data: gadeduc.csv

This is synthetic data from a randomised controlled trial, in which 30 therapists randomly assigned patients (each therapist saw between 2 and 28 patients) to a control or treatment group, and monitored their scores over time on a measure of generalised anxiety disorder (GAD7 - a 7 item questionnaire with 5 point likert scales).

The control group of patients received standard sessions offered by the therapists. For the treatment group, 10 mins of each sessions was replaced with a specific psychoeducational component, and patients were given relevant tasks to complete between each session. All patients had monthly therapy sessions. Generalised Anxiety Disorder was assessed at baseline and then every visit over 4 months of sessions (5 assessments in total).

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

You can find a data dictionary below:

Table 1: Data Dictionary: lmm_gadeduc.csv

variable	description
patient	A patient code in which the labels take the form <Therapist initials>_<group>_<patient number>.
visit_0	Score on the GAD7 at baseline
visit_1	GAD7 at 1 month assessment
visit_2	GAD7 at 2 month assessment
visit_3	GAD7 at 3 month assessment
visit_4	GAD7 at 4 month assessment

Question 1

Uh-oh… these data aren’t in the same shape as the other datasets we’ve been giving you..

Can you get it into a format that is ready for modelling?

Hints

• It’s wide, and we want it long.
  
• Once it’s long. “visit_0”, “visit_1”,.. needs to become the numbers 0, 1, …
• One variable (patient) contains lots of information that we want to separate out. There’s a handy function in the tidyverse called separate(), check out the help docs!

Solution 1. Here’s the data. We have one row per patient, but we have multiple observations for each patient across the columns..

geduc = read_csv("https://uoepsy.github.io/data/lmm_gadeduc.csv")
head(geduc)

# A tibble: 6 × 6
  patient      visit_0 visit_1 visit_2 visit_3 visit_4
  <chr>          <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
1 VC_Control_1      24      24      26      29      28
2 VC_Control_2      24      26      28      29      30
3 VC_Control_3      25      29      27      29      30
4 VC_Control_4      24      25      25      26      26
5 VC_Control_5      28      28      27      29      28
6 VC_Control_6      26      28      25      27      28

We can make it long by taking the all the columns from visit_0 to visit_4 and shoving their values into one variable, and keeping the name of the column they come from as another variable:

geduc |> 
  pivot_longer(2:last_col(), names_to="visit",values_to="GAD")

# A tibble: 2,410 × 3
   patient      visit     GAD
   <chr>        <chr>   <dbl>
 1 VC_Control_1 visit_0    24
 2 VC_Control_1 visit_1    24
 3 VC_Control_1 visit_2    26
 4 VC_Control_1 visit_3    29
 5 VC_Control_1 visit_4    28
 6 VC_Control_2 visit_0    24
 7 VC_Control_2 visit_1    26
 8 VC_Control_2 visit_2    28
 9 VC_Control_2 visit_3    29
10 VC_Control_2 visit_4    30
# ℹ 2,400 more rows

Solution 2. Now we know how to get our data long, we need to sort out our time variable (visit) and make it into numbers.
We can replace all occurrences of the string "visit_" with nothingness "", and then convert them to numeric.

geduc |> 
  pivot_longer(2:last_col(), names_to="visit",values_to="GAD") |>
  mutate(
    visit = as.numeric(gsub("visit_","",visit))
  ) 

# A tibble: 2,410 × 3
   patient      visit   GAD
   <chr>        <dbl> <dbl>
 1 VC_Control_1     0    24
 2 VC_Control_1     1    24
 3 VC_Control_1     2    26
 4 VC_Control_1     3    29
 5 VC_Control_1     4    28
 6 VC_Control_2     0    24
 7 VC_Control_2     1    26
 8 VC_Control_2     2    28
 9 VC_Control_2     3    29
10 VC_Control_2     4    30
# ℹ 2,400 more rows

Solution 3. Finally, we need to sort out the patient variable. It contains 3 bits of information that we will want to have separated out. It has the therapist (their initials), then the group (treatment or control), and then the patient number. These are all separated by an underscore “_“.

The separate() function takes a column and separates it into several things (as many things as we give it), splitting them by some user defined separator such as an underscore:

geduc_long <- geduc |> 
  pivot_longer(2:last_col(), names_to="visit",values_to="GAD") |>
  mutate(
    visit = as.numeric(gsub("visit_","",visit))
  ) |>
  separate(patient, into=c("therapist","group","patient"), sep="_")

And we’re ready to go!

geduc_long

# A tibble: 2,410 × 5
   therapist group   patient visit   GAD
   <chr>     <chr>   <chr>   <dbl> <dbl>
 1 VC        Control 1           0    24
 2 VC        Control 1           1    24
 3 VC        Control 1           2    26
 4 VC        Control 1           3    29
 5 VC        Control 1           4    28
 6 VC        Control 2           0    24
 7 VC        Control 2           1    26
 8 VC        Control 2           2    28
 9 VC        Control 2           3    29
10 VC        Control 2           4    30
# ℹ 2,400 more rows

Question 2

Visualise the data. Does it look like the treatment had an effect?
Does it look like it worked for every therapist?

Hints

• remember, stat_summary() is very useful for aggregating data inside a plot.

Here’s the overall picture. The average score on the GAD7 at each visit gets more and more different between the two groups. The treatment looks effective..

ggplot(geduc_long, aes(x = visit, y = GAD, col = group)) +
  stat_summary(geom="pointrange")

Let’s split this up by therapist, so we can see the averages across each therapist’s set of patients.
There’s clear variability between therapists in how well the treatment worked. For instance, the therapists EU and OD don’t seem to have much difference between their groups of patients.

ggplot(geduc_long, aes(x = visit, y = GAD, col = group)) +
  stat_summary(geom="pointrange") +
  facet_wrap(~therapist)

Question 3

Fit a model to test if the psychoeducational treatment is associated with greater improvement in anxiety over time.

We want to know if how anxiety (GAD) changes over time (visit) is different between treatment and control (group).

Hopefully this should hopefully come as no surprise¹ - it’s an interaction!

lmer(GAD ~ visit * group + ...
       ...
     data = geduc_long)

We have multiple observations for each of the 482 patients, and those patients are nested within 30 therapists.

Note that in our data, the patient variable does not uniquely specify the individual patients. i.e. patient “1” from therapist “AO” is a different person from patient “1” from therapist “BJ”. To correctly group the observations into different patients (and not ‘patient numbers’), we need to have therapist:patient.

So we capture therapist-level differences in ( ... | therapist) and the patients-within-therapist-level differences in ( ... | therapist:patient):

lmer(GAD ~ visit * group + ...
       ( ... | therapist) + 
       ( ... | therapist:patient),
     data = geduc_long)

Note that each patient can change differently in their anxiety levels over time - i.e. the slope of visit could vary by participant.

Likewise, some therapists could have patients who change differently from patients from another therapist, so visit|therapist can be included.

Each patient is in one of the two groups - they’re either treatment or control. So we can’t say that “differences in anxiety due to treatment varies between patients”, because for any one patient the “difference in anxiety due to treatment” is not defined in our study design.

However, therapists see multiple different patients, some of which are in the treatment group, and some of which are in the control group. So the treatment effect could be different for different therapists!

mod1 <- lmer(GAD ~ visit*group + 
               (1+visit*group|therapist)+
               (1+visit|therapist:patient),
             geduc_long)

Question 4

For each of the models below, what is wrong with the random effect structure?

modelA <- lmer(GAD ~ visit*group + 
               (1+visit*group|therapist)+
               (1+visit|patient),
             geduc_long)

modelB <- lmer(GAD ~ visit*group + 
               (1+visit*group|therapist/patient),
             geduc_long)

modelA <- lmer(GAD ~ visit*group + 
               (1+visit*group|therapist)+
               (1+visit|patient),
             geduc_long)

The patient variable doesn’t capture the different patients within therapists, so this actually fits crossed random effects and treats all data where patient==1 as from the same group (even if this includes several different patients’ worth of data from different therapists!)

modelB <- lmer(GAD ~ visit*group + 
               (1+visit*group|therapist/patient),
             geduc_long)

Using the / here means we have the same random slopes fitted for therapists and for patients-within-therapists. but the effect of group can’t vary by patient, so this doesn’t work. hence why we need to split them up into (...|therapist)+(...|therapist:patient).

Question 5

Let’s suppose that I don’t want the psychoeducation treatment, I just want the standard therapy sessions that the ‘Control’ group received. Which therapist should I go to?

Hints

dotplot.ranef.mer() might help here!
You can read about ranef in Chapter 2 #making-model-predictions.

It would be best to go to one of the therapists SZ, YS, or IT…

Why? These therapists all have the most negative slope of visit:

dotplot.ranef.mer(ranef(mod1))$therapist

Question 6

Recreate this plot.

The faint lines represent the model estimated lines for each patient. The points and ranges represent our fixed effect estimates and their uncertainty.

Hints

• you can get the patient-specific lines using augment() from the broom.mixed package, and the fixed effects estimates using the effects package.
• remember that the “patient” column doesn’t group observations into unique patients.
• remember you can pull multiple datasets into ggplot:

ggplot(data = dataset1, aes(x=x,y=y)) + 
  geom_point() + # points from dataset1
  geom_line(data = dataset2) # lines from dataset2

• see more in Chapter 2 #visualising-models

Solution 4. The effects package will give us the fixed effect estimates:

library(effects)
library(broom.mixed)
effplot <- effect("visit*group",mod1) |>
  as.data.frame()

We want to get the fitted values for each patient. We can get fitted values using augment(). But the patient variable doesn’t capture the unique patients, it just captures their numbers (which aren’t unique to each therapist).
So we can create a new column called upatient which pastes together the therapists initials and the patient numbers

augment(mod1) |> 
  mutate(
    upatient = paste0(therapist,patient),
    .after = patient # place the column next to the patient col
  )

# A tibble: 2,410 × 17
     GAD visit group   therapist patient upatient .fitted .resid  .hat .cooksd
   <dbl> <dbl> <fct>   <fct>     <fct>   <chr>      <dbl>  <dbl> <dbl>   <dbl>
 1    24     0 Control VC        1       VC1         24.2 -0.198 0.454 0.0210 
 2    24     1 Control VC        1       VC1         25.3 -1.28  0.239 0.239  
 3    26     2 Control VC        1       VC1         26.4 -0.360 0.186 0.0128 
 4    29     3 Control VC        1       VC1         27.4  1.56  0.294 0.508  
 5    28     4 Control VC        1       VC1         28.5 -0.522 0.563 0.284  
 6    24     0 Control VC        2       VC2         24.8 -0.843 0.454 0.383  
 7    26     1 Control VC        2       VC2         26.2 -0.171 0.239 0.00426
 8    28     2 Control VC        2       VC2         27.5  0.502 0.186 0.0250 
 9    29     3 Control VC        2       VC2         28.8  0.174 0.294 0.00633
10    30     4 Control VC        2       VC2         30.2 -0.153 0.563 0.0246 
# ℹ 2,400 more rows
# ℹ 7 more variables: .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>,
#   .sqrtrwt <dbl>, .weights <dbl>, .wtres <dbl>

Solution 5.

library(effects)
library(broom.mixed)
effplot <- effect("visit*group",mod1) |>
  as.data.frame()

augment(mod1) |> 
  mutate(
    upatient = paste0(therapist,patient),
    .after = patient # place the column next to the patient col
  ) |>
  ggplot(aes(x=visit,y=.fitted,col=group))+
  stat_summary(geom="line", aes(group=upatient,col=group), alpha=.1)+
  geom_pointrange(data=effplot, aes(y=fit,ymin=lower,ymax=upper,col=group))+
  labs(x="- Month -",y="GAD7")

Jokes

Data: lmm_laughs.csv

These data are simulated to imitate an experiment that investigates the effect of visual non-verbal communication (i.e. gestures, facial expressions) on joke appreciation. 90 Participants took part in the experiment, in which they each rated how funny they found a set of 30 jokes. For each participant, the order of these 30 jokes was randomly set for each run of the experiment. For each participant, the set of jokes was randomly split into two halves, with the first half being presented in audio-only, and the second half being presented in audio and video. This meant that each participant saw 15 jokes with video and 15 without, and each joke would be presented in with video roughly half of the times it was seen.

The researchers want to investigate whether the delivery (audio/audiovideo) of jokes is associated with differences in humour-ratings.

Data are available at https://uoepsy.github.io/data/lmm_laughs.csv

Table 2: Data Dictionary: lmm_laughs.csv

variable	description
ppt	Participant Identification Number
joke_label	Joke presented
joke_id	Joke Identification Number
delivery	Experimental manipulation: whether joke was presented in audio-only ('audio') or in audiovideo ('video')
rating	Humour rating chosen on a slider from 0 to 100

Question 7

Prior to getting hold of any data or anything, we should be able to write out the structure of our ideal “maximal” model.

Can you do so?

Hints

Don’t know where to start? Try following the steps in Chapter 8 #maximal-model.

We want to estimate the effect of delivery on humour rating of jokes:
rating ~ delivery

We have 30 observations for each participant. Participants are just another sampling unit here.
rating ~ delivery + (1 | ppt)

We have 90 observations for each joke. We’re not interested in specific jokes here, so we can think of these as a random set of experimental items that we might choose differently next time we conduct an experiment to assess delivery~rating: rating ~ delivery + (1 | ppt) + (1 | joke_id)

Participants each see 15 jokes without video, and 15 with. The delivery variable is “within” participant. Some participants might respond a lot to having video whereas some might not rate jokes any differently. The effect of delivery on rating might be vary by participant:
rating ~ delivery + (1 + delivery | ppt) + (1 | joke_id)

Each joke is presented both with and without the video. Some jokes might really benefit from gestures and facial expressions, whereas some might not. The effect of delivery on rating might be vary by joke:
rating ~ delivery + (1 + delivery | ppt) + (1 + delivery | joke_id)

Question 8

Read in and clean the data (if necessary).

Create some plots showing:

1. the average rating for audio vs audio+video for each joke
2. the average rating for audio vs audio+video for each participant

Hints

• you could use facet_wrap, or even stat_summary!
  
• you might want to use joke_id, rather than joke_label (the labels are very long!)

Here is one using facet_wrap:

ggplot(laughs, aes(x = delivery, y = rating)) +
  geom_boxplot()+
  facet_wrap(~joke_id)

And one using stat_summary() for participants:

ggplot(laughs, aes(x = delivery, y = rating)) +
  stat_summary(geom="pointrange", aes(group = ppt),
               position = position_dodge(width=.2))+
  stat_summary(geom="line", aes(group = ppt),
               position = position_dodge(width=.2))

Question 9

Fit an appropriate model to address the research aims of the study.

Hints

This should be the one you came up with a couple of questions ago!

mod <- lmer(rating ~ delivery + 
              (1 + delivery | joke_id) +
              (1 + delivery| ppt), data = laughs)

Question 10

Which joke is funniest when presented just in audio? For which joke does the video make the most difference to ratings?

Hints

These can all be answered by examining the random effects with ranef().
See Chapter 2 #making-model-predictions.

If you’re using joke_id, can you find out the actual joke that these correspond to?

dotplot.ranef.mer(ranef(mod))$joke_id

Joke 19 is the funniest apparently! (not sure I agree)

Lots of ways to find what the joke actually is. Here is one way:

laughs |> count(joke_id, joke_label) |>
  filter(joke_id==19) |>
  pull(joke_label)

[1] "How many psychiatrists does it take to change a lightbulb? Just one, but the lightbulb really has to want to change."

And from the plot above, Joke 28 has the most benefit of video. We can quickly check this with something like:

ranef(mod)$joke_id |>
  filter(deliveryvideo == max(deliveryvideo))

   (Intercept) deliveryvideo
28       0.807          3.37

The joke itself is a bit weird.. maybe the video really helped!

laughs |> count(joke_id, joke_label) |>
  filter(joke_id==28) |>
  pull(joke_label)

[1] "An Alsatian went to a telegram office, took out a blank form and wrote:\n\"Woof. Woof. Woof. Woof. Woof. Woof. Woof. Woof. Woof.\"\nThe clerk examined the paper and politely told the dog: \"There are only nine\nwords here. You could send another \x91Woof' for the same price.\"\n\"But,\" the dog replied, \"that would make no sense at all.\""

Question 11

Do jokes that are rated funnier when presented in audio-only tend to also benefit more from the addition of video?

Hints

Think careful about this question. The random effects show us that jokes vary in their intercepts (ratings in audio-only) and in their effects of delivery (the random slopes). We want to know if these are related..

VarCorr(mod)

 Groups   Name          Std.Dev. Corr
 ppt      (Intercept)   3.51         
          deliveryvideo 1.78     0.00
 joke_id  (Intercept)   2.15         
          deliveryvideo 2.24     0.39
 Residual               5.83         

It’s the correlation here that tell us - jokes rated higher in the audio-only tend to have a bigger effect of the video.

We can see this in a plot if we like. Here every dot is a joke, and the x-axis shows whether it is above or below the average rating for audio-only (the intercept). The y-axis shows whether it is above or below the average effect of video.

plot(ranef(mod)$joke)

Question 12

Create a plot of the estimated effect of video on humour-ratings. Try to plot not only the fixed effects, but the raw data too.

Hints

See e.g. Chapter 2 #visualising-models

library(effects)

plotdatf <- effect("delivery",mod) |>
  as.data.frame()

ggplot(data = laughs, aes(x = delivery)) +
  geom_jitter(aes(y = rating), width = .1, height = 0, alpha = .1) +
  geom_pointrange(data = plotdatf,
                  aes(y = fit, ymin = lower, ymax = upper),
                  position=position_nudge(x=.2))

Extra: Test Enhanced Learning

Data: Test-enhanced learning

An experiment was run to conceptually replicate “test-enhanced learning” (Roediger & Karpicke, 2006): two groups of 25 participants were presented with material to learn. One group studied the material twice (StudyStudy), the other group studied the material once then did a test (StudyTest). Recall was tested immediately (one minute) after the learning session and one week later. The recall tests were composed of 175 items identified by a keyword (Test_word).

The critical (replication) prediction is that the StudyStudy group perform better on the immediate test, but the StudyTest group will retain the material better and thus perform better on the 1-week follow-up test.

Test performance is measured as the speed taken to correctly recall a given word.

The following code loads the data into your R environment by creating a variable called tel:

load(url("https://uoepsy.github.io/data/testenhancedlearning.RData"))

Table 3: Data Dictionary: testenhancedlearning.Rdata

variable	description
Subject_ID	Unique Participant Identifier
Group	Group denoting whether the participant studied the material twice (StudyStudy), or studied it once then did a test (StudyTest)
Delay	Time of recall test ('min' = Immediate, 'week' = One week later)
Test_word	Word being recalled (175 different test words)
Correct	Whether or not the word was correctly recalled
Rtime	Time to recall word (milliseconds)

Question 13

Here is the beginning of our modelling.

# load in the data
load(url("https://uoepsy.github.io/data/testenhancedlearning.RData"))

# performance is measured by time taken to *correctly*
# recall a word.
# so we're going to have to discard all the incorrects:
tel <- tel |> filter(Correct == 1)

# preliminary plot
# makes it look like studytest are better at immediate (contrary to prediction)
# both groups get slower from immediate > week, 
# but studystudy slows more than studytest
ggplot(tel, aes(x = Delay, y = Rtime, col = Group)) + 
  stat_summary(geom="pointrange")

mmod <- lmer(Rtime ~ Delay*Group +
             (1 + Delay | Subject_ID) +
             (1 + Delay * Group | Test_word),
             data=tel)

This is what I did. You might do something else!

First I removed the interaction from the random effects

mod1 <- lmer(Rtime ~ Delay*Group +
             (1 + Delay | Subject_ID) +
             (1 + Delay + Group | Test_word),
             data=tel)

boundary (singular) fit: see help('isSingular')

This model is a singular fit, suggesting it needs further simplification. The variance covariance matrix of the random effects isn’t giving us many pointers..

# examine vcov
VarCorr(mod1)

 Groups     Name           Std.Dev. Corr       
 Test_word  (Intercept)     21.5               
            Delayweek       14.5     0.17      
            GroupStudyTest  23.4    -0.77 -0.76
 Subject_ID (Intercept)     27.1               
            Delayweek       10.2    0.01       
 Residual                  240.1               

There are various things we could try here. See Chapter 8 #simplifying-random-effect-structures for some of the more in-depth options.

However, sometimes it is simplest just to trial & error the removal of different possible terms. Here we are removing Delay|Test_word and removing Delay|Subject_ID:

mod2a <- lmer(Rtime ~ Delay*Group +
             (1 + Delay | Subject_ID) +
             (1 + Group | Test_word),
             data=tel)
mod2b <- lmer(Rtime ~ Delay*Group +
             (1 | Subject_ID) +
             (1 + Delay + Group | Test_word),
             data=tel)

boundary (singular) fit: see help('isSingular')

The second model is a singular fit, but the first one is not. Just for safety, let’s check:

isSingular(mod2a)

[1] FALSE

All looks good there.

Sometimes it can be useful to check how estimates of fixed effects and their standard errors differ across possible candidate models with different random effect structures. More often than not, this simply provides us with reassurance that the removal of random effects hasn’t actually had too much of an impact on anything we’re going to conduct inferences on. If they differ a lot, then we have a lot more to discuss!

Here are the fixed effects from each model:

term	mod1	mod2a	mod2b
(Intercept)	740.55 (7.17)	740.57 (7.21)	740.69 (7.23)
Delayweek	27.65 (6.97)	27.64 (6.87)	27.23 (6.64)
GroupStudyTest	-31.82 (10.26)	-31.75 (10.23)	-31.73 (10.35)
Delayweek:GroupStudyTest	-17.2 (9.69)	-17.26 (9.7)	-17.18 (9.19)

In all these models, the fixed effects estimates are all pretty similar, suggesting that they’ve all found similar estimates of these parameters which have been largely invariant to our refinement of the random effects. This makes me feel better - there’s less worry that our final conclusions are going to be influenced by specifics of incl/exclusion of one of these random effect terms.

I would definitely settle on mod2a because that is the one that converges, but we can add a footnote if we wanted, to mention that mod2b finds the same pattern of results.

Footnotes

1. if it does, head back to where we learned about interactions in the single level regressions lm(). It’s just the same here.