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 https://stats.stackexchange.com/questions/639430/frank-harrells-interpretation-of-interaction-in-regression-results/639444#639444. 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.
```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"
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"
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. 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 and should not be done by the regular likelihood-ratio test such as through `anova()`. It requires specific likelihood-ratio test methods for variance components, usually included in the package.