(Mixed-)model specification and comparison for testing hypotheses concerning continuous ratings

Question

In a behavioral study following a within-subjects design, subjects of various levels of musicality (years of education) gave continuous liking ratings of three audio stimuli, each under three listening conditions: M, L, and ML. Ratings were recorded continuously, with one sample encoded from a rating slider every 0.5s.

The data is stored in a table where every row refers to one sample, with the following explanatory variables to be defined as fixed or random effects:

• [DV: Rating given]
• Stimulus number
• Condition
• Subject musicality
• Subject number
• Sample number
• Loudness

The research question concerns the difference between average ratings under the three conditions. Specific hypotheses to be tested are:

1. ML ratings contain information that is independent of any straightforward "combination" of M ratings and L ratings; this might be more the case for more musical subjects

2. ML ratings will more closely be followed by M ratings than by L ratings; again, this might hold more for more musical subjects

3. The three stimuli will not differ qualitatively from one another in terms of the degree to which hypotheses 1 and 2 hold 

Being new to LMMs, I am unsure how to proceed, namely:

• how to define the various linear mixed models (including any null/base models if needed) and their random intercepts and/or slopes
• which of these models to compare* in order to test each hypothesis (* and whether to base said comparisons on likelihood ratio tests or on AICs)
• how to account for the serial correlation within each rating signal, i.e. the fact that any one sample won't be massively different from the previous sample or the next one, given that ratings are made using a physical slider

I realise this is a complex question but I'm not sure how else to simplify it. Any code to implement this in Matlab would be greatly appreciated as an extra!

Answer 1

According to your comment that "One measurement every 0.5s, and this was done for 3 separate stimuli, each under 3 conditions, so 9 in total," I assume that each Stimulus corresponded to a unique song, and three songs were selected in this experiment. Whereas each Subject number identifies a person, each Sample number is an integer serial stamping each half second elapsed from the beginning of a song. The frequency of two ratings per second sounds really weird and unrealistic though. I doubt if subjects had enough time to evaluate to provide a meaningful response. If each song lasted three minutes, this means that each Subject gave 3 * 120 ratings to each Stimulus × Condition combination, so there were 3240 measurements from each person. High-frequency data, such as stock trading series, may not be sufficiently represented in mixed effects models. You may need time series analysis instead.

If the ratings were limited categorical levels, such as 1 = very boring, 5 = very interesting, I suggest trying ordinal regression with random effects. See R package {ordinal} https://cran.r-project.org/web/packages/ordinal/index.html and specifically Cumulative Link Mixed Model https://cran.r-project.org/web/packages/ordinal/vignettes/clmm2_tutorial.pdf.

I suggest relabeling the sample index as Time or Timestamp because a sample usually refers to a collection of measurements in statistics. It may also be beneficial to treat this time variable continuous instead of an index. You could build Fourier terms to see how ratings change over a song's playtime. If monotonic, you could use exponentiation or logarithm to represent time decay, such as Rating ~ log(Time).

According to your hypothesis that ML is close to M and higher than L, a convenient encoding should consider L as the reference group. I suggest breaking the Condition variable into two indicators, LTrue and MTrue, to represent whether the song was played under M and whether the song was played under L, respectively. See the meaning of coefficients in a common model specification with an interaction Frank Harrell's interpretation of interaction in regression results. Usually two binary indicators form four groups, but here you had only three, missing the group that song was played under neither M nor L. With proper design matrix manipulation via model formula, we can easily see how much higher effect ML brings than M, which to be contrasted with L. To remove potential redundant terms which would otherwise cause perfect multicollinearity, we may need to build specific interaction terms using : instead of * in R. See a comparison of different ways to specifying the same model in R.

library(dplyr)
Data <- data.frame(
  Rating = 1:30, 
  Condition = factor(c(rep("L", 10), rep("M", 10), rep("LM", 10)), levels = c(
    "L", "M", "LM")), 
  L01 = c(rep(1, 10), rep(0, 10), rep(1, 10)),
  M01 = c(rep(0, 10), rep(1, 10), rep(1, 10)), 
  LTrue = c(rep(TRUE, 10), rep(FALSE, 10), rep(TRUE, 10)),
  MTrue = c(rep(FALSE, 10), rep(TRUE, 10), rep(TRUE, 10)), 
  Lny = factor(c(rep("Yes", 10), rep("No", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")),
  Mny = factor(c(rep("No", 10), rep("Yes", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")))
Data |>
  group_by(Condition) |>
  summarise(mean = mean(Rating))
"# A tibble: 3 × 2
  Condition  mean
  <fct>     <dbl>
1 L           5.5
2 M          15.5
3 LM         25.5"

# Factor treatment contrast
summary(lm(Rating ~ Condition, data = Data))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
ConditionM   10.0000     1.3540   7.385 6.05e-08 ***
ConditionLM  20.0000     1.3540  14.771 1.87e-14 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# 0/1 encoding
summary(lm(Rating ~ 1 + L01 * M01, data = Data)) # weird coef
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.500      1.658  -2.714   0.0114 *  
L01           10.000      1.354   7.385 6.05e-08 ***
M01           20.000      1.354  14.771 1.87e-14 ***
L01:M01           NA         NA      NA       NA
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + L01 * M01, data = Data)) # wrong R2 
"        Estimate Std. Error t value Pr(>|t|)    
L01       5.5000     0.9574   5.745 4.16e-06 ***
M01      15.5000     0.9574  16.189 2.00e-15 ***
L01:M01   4.5000     1.6583   2.714   0.0114 *
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ M01 + L01 : M01, data = Data)) # correct formula
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
M01          10.0000     1.3540   7.385 6.05e-08 ***
M01:L01      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Logical encoding
summary(lm(Rating ~ 1 + LTrue * MTrue, data = Data)) # weird coef
"                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -4.500      1.658  -2.714   0.0114 *  
LTrueTRUE             10.000      1.354   7.385 6.05e-08 ***
MTrueTRUE             20.000      1.354  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + LTrue * MTrue, data = Data)) # weird coef + wrong R2
"                    Estimate Std. Error t value Pr(>|t|)    
LTrueFALSE           -4.5000     1.6583  -2.714   0.0114 *  
LTrueTRUE             5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE            20.0000     1.3540  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ MTrue + LTrue : MTrue, data = Data)) # correct formula
"                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE             10.0000     1.3540   7.385 6.05e-08 ***
MTrueFALSE:LTrueTRUE       NA         NA      NA       NA    
MTrueTRUE:LTrueTRUE   10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ MTrue + I(LTrue * MTrue), data = Data)) # Removes NA above
"                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE         10.0000     1.3540   7.385 6.05e-08 ***
I(LTrue * MTrue)  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Binary factor encoding
summary(lm(Rating ~ 1 + Lny * Mny, data = Data)) # weird coef
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -4.500      1.658  -2.714   0.0114 *  
LnyYes          10.000      1.354   7.385 6.05e-08 ***
MnyYes          20.000      1.354  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + Lny * Mny, data = Data)) # weird coef + wrong R2
"              Estimate Std. Error t value Pr(>|t|)    
LnyNo          -4.5000     1.6583  -2.714   0.0114 *  
LnyYes          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         20.0000     1.3540  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ Mny + Lny : Mny, data = Data)) # correct formula
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         10.0000     1.3540   7.385 6.05e-08 ***
MnyNo:LnyYes        NA         NA      NA       NA    
MnyYes:LnyYes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny : Mny), data = Data)) # NA still there
"                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.5000     0.9574   5.745 4.16e-06 ***
MnyYes             10.0000     1.3540   7.385 6.05e-08 ***
I(Lny:Mny)Yes:No        NA         NA      NA       NA    
I(Lny:Mny)Yes:Yes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny == "Yes" & Mny == "Yes"), data = Data)) # no NA
"                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes                              10.0000     1.3540   7.385 6.05e-08 ***
I(Lny == Yes & Mny == Yes)TRUE      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

Notice that using a three-level factor gives coefficients as contrasts to L. Separating the factor into two binary indicators, however, often introduces in redundant terms if there is a missing group. Removing the intercept to combat multicollinearity often jeopardizes some model summary statistics like $R^2$. Instead, specifying one main effect and one interaction M + L : M seem to offer the easiest interpretation among three groups where the coefficient of M gives its higher effect than L and the coefficient of L : M interaction gives the difference between ML and M. To find the difference between ML and L + M, use the difference between the coefficient of L : M interaction and the intercept.

which of these models to compare* in order to test each hypothesis (* and whether to base said comparisons on likelihood ratio tests or on AICs)

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American Statistician, 61(1), 22–27. https://www.jstor.org/stable/27643833.

Answer 2

To address the research questions and hypotheses using a mixed effects approach modelling approach, we need to carefully specify the fixed and random effects. Based on your research question, the obvious fixed effects are:

• Stimulus number (stimulus in what follows below)
• Condition (condition in what follows below)
• Subject musicality (musicality)
• Sample number (sample)
• Loudness (loudness)

The obvious candidate for random intercepts is the subject. You don't say how many subjects you have. There is no hard and fast rule for this but general consensus (in my experience says) that fewer than 6 subjects may be problematic.

A baseline model could be:

rating ~ stimulus + condition + musicality + sample + loudness + (1|subject)

That is, main effects (no interactions) for the fixed effects and random intercepts for subjects. Your research questions can (mostly) be addressed by extending this model with interactions:

rating ~ stimulus * condition * musicality + sample + loudness + (1|subject)

1. ML ratings contain information that is independent of any straightforward "combination" of M ratings and L ratings; this might be more the case for more musical subjects.

If the two-way interaction between condition and musicality is significant, it indicates that the effect of the condition (M, L, ML) on the rating is modulated by the subject's musicality.

If the interaction between the ML condition and musicality is significant and distinct from the sum of coefficients for the M (if M is used as the baseline or reference condition, its coefficient ought to be zero) and L conditions, it supports the hypothesis that ML ratings contain unique information that is not just a combination of M and L ratings. This effect being stronger for more musical subjects would be reflected in a larger coefficient of the interaction term between ML and musicality than that of the interaction term between L and musicality.

2. ML ratings will more closely be followed by M ratings than by L ratings; again, this might hold more for more musical subjects.

This hypothesis can be partially assessed by examining how the ratings under the ML condition relate to subsequent ratings under the M condition, especially in terms of interaction with musicality. While the model does not directly model time-lagged effects, the significance and direction of the condition:interaction terms can provide insight. If the ML condition's effect is significantly different from both M and L conditions, and the interaction with musicality is significant, it suggests that ML ratings are uniquely related to M ratings for more musical subjects. Obviously this does not fully answer the question. To fully address this question, you may need to extend the model to explicitly include lagged effects, and we can do so as follows:

rating ~ rating_prev*condition*musicality + stimulus*condition*musicality + sample + loudness + (1|subject)

This answers the research question by including interactions between the lagged rating and the condition and musicality. This allows us to determine if the previous rating (under different conditions) influences the current rating more strongly when the previous condition was ML compared to M or L, and whether this effect is modulated by the subject's musicality.

Obviously you will need to modify your dataframe to incorporate the lagged ratings.

3. The three stimuli will not differ qualitatively from one another in terms of the degree to which hypotheses 1 and 2 hold.

If the three-way interaction term stimulus:condition:musicality is not significant, it suggests that the stimuli do not differ regarding the effects described in questions 1 and 2.

You also ask:

how to account for the serial correlation within each rating signal, i.e. the fact that any one sample won't be massively different from the previous sample or the next one, given that ratings are made using a physical slider

One way to account for serial correlation is to fit random slopes for the time variable:

rating ~ stimulus * condition * musicality + sample + loudness + (sample|subject)

Another way to accomplish it is to specify an autoregressive residual error covariance structure. As far as I know MATLAB's mixed model procedure does not provide this functionality. In R, you could use glmmTMB, nlme or the mmrm packages.

I should mention that the model I suggested above with the lagged rating also handles the autoregression. You might want to consider using that model to answer ALL your research questions. I hesitate to suggest this because the model was already quite complex, with a 3-way interaction; the lagged model has two 3-way interactions, and then when you include random slopes there will be an avalanche of output to interpret - and that assumes that such a complex model will even converge normally. If it does, and you are happy working with the lagged model to answer all the questions, then by all means do that !

And, one further point - since the measurements are so close to each other (0.5s), a compound symmetric residual covariance structure might be fine - this means that the correlation between any two observations in time are the same - regardless of the time between them. This structure is induced by a simple random intercepts model, so you might find that using random slopes or an AR model might not be needed.

Any code to implement this in Matlab would be greatly appreciated as an extra!

My MATLAB skills are rather rusty, but hopefully, this will get you going:

% Load your data
data = readtable('your_data.csv');

% Base Model
baseModel = fitlme(data, 'rating ~ stimulus + condition + musicality + sample + loudness + (1|subject)');

% Interaction Effects Model
interactionModel = fitlme(data, 'rating ~ stimulus * condition * musicality + sample + loudness + (1 |subject)');

% Random Slopes Model
randomSlopesModel = fitlme(data, 'rating ~ stimulus * condition * musicality + sample + loudness + (1 + sample|subject)');

% Model Comparison using AIC
AIC_values = [baseModel.ModelCriterion.AIC, interactionModel.ModelCriterion.AIC, randomSlopesModel.ModelCriterion.AIC];
disp('AIC values for models:');
disp(AIC_values);

% Model Comparison using Likelihood Ratio Test
compare(interactionModel, baseModel)
compare(randomSlopesModel, interactionModel)

One final note: You don't provide much detail about the response variable - above I have assumed that it is continuous. It will be important to check that the residuals are plausibly normally distributed - if not, then you might need to transform it, or use an appropriate generalized linear mixed model (GLMM).