# Mixed effects model to explore the linear relationship between two variables for each level of the repeated-measures in Rstudio using lme()

Background: I did an experiment in which 30 people watched 3 movie clips ("A", "B", and "C") in a random order. For each song and participant I measured one continuous neural variable named SNR (1 value for each participant and for each song) and one continuous behavioral variable named RT (1 value for each participant and for each song). I am interested in knowing whether there is a linear relationship between SNR and RT for each movie. Also, I am interested in knowing the strength of this relationship to compare across movies. In another words, I want to know whether the greater the SNR the greater the RT for each movie clip, and whether this is enhanced in some movies compared to others.

Here is some dummy data that follows the same structure that I have:

res_table <- tibble(
ID                   = rep(seq(1,30,1),3), # participant identifier
MOVIES               = rep(c("A","B","C"),each=30),
SNR                  = sample(35:65, 30*3, replace = TRUE),
RT                   = sample(50:100, 30*3, replace = TRUE)
)
res_table$$ID <- factor(res_table$$ID, labels = c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"))
res_table$$MOVIES <- factor(res_table$$MOVIES, levels = c("A","B","C"))


Problem: Although I initially performed Pearson Correlations between SNR and RT for each movie, someone told me that to explore these relationships it would be best to "show a scatter plot of the association between SNR and RT, including the linear model and the confidence range for the regression".

I fitted a mixed model as follows:

model.01 <- lmer(SNR~RT*MOVIES + (1|ID), data = res_table)


But from here on I am lost on what I should do. How can I know whether the relationship between SNR and RT at each movie level is significant or not?

Using the data you provided: look at the model summary first:

summary(model.01)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: SNR ~ RT * MOVIES + (1 | ID)
Data: res_table

REML criterion at convergence: 636

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.5020 -0.7392 -0.1719  0.8393  1.7881

Random effects:
Groups   Name        Variance Std.Dev.
ID       (Intercept) 13.03    3.610
Residual             62.24    7.889
Number of obs: 90, groups:  ID, 30

Fixed effects:
Estimate Std. Error       df t value Pr(>|t|)
(Intercept) 38.51328    6.96568 80.81879   5.529 3.84e-07 ***
RT           0.15065    0.09162 80.03000   1.644  0.10404
MOVIESB     29.10257   10.32089 77.10354   2.820  0.00611 **
MOVIESC     16.89379   10.63170 80.86384   1.589  0.11596
RT:MOVIESB  -0.38552    0.13743 77.85638  -2.805  0.00635 **
RT:MOVIESC  -0.24316    0.13840 81.47951  -1.757  0.08268 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The summary suggests that the relation between RT and SNR is different for movie A vs. movie B at least. However you probably want to compare the slope of RT for each different movie. You can do this via emmeans:

library(emmeans)
em<-emtrends(model.01, "MOVIES", var="RT")
em
MOVIES RT.trend     SE   df lower.CL upper.CL   #RT.trend gives you the RT coefficient
A        0.1507 0.0934 80.2  -0.0352   0.3365  #separately for each movie
B       -0.2349 0.1058 80.3  -0.4454  -0.0244
C       -0.0925 0.1047 80.3  -0.3008   0.1158

#then you can compare the trends by

pairs(em)
contrast estimate    SE   df t.ratio p.value
A - B       0.386 0.140 78.2   2.756  0.0197
A - C       0.243 0.141 81.6   1.722  0.2033
B - C      -0.142 0.151 82.5  -0.945  0.6135


This suggests that 1) RT is only predicting SNR significantly for movie type B, and 2) RT effect on SNR for movie A is statistically significantly different from RT effect on SNR for movie B (at p < .05), however there is no significant difference in RT effect between A and C or B and C.