# how to specify a hierarchical mixed model with two factors nested within participant that are crossed with each other

I've read a bunch of questions and answers here about how to specify models and some tutorial walk throughs and I still not sure how to deal with my data structure. I want, obviously, to avoid any hint of pseudoreplication.

I'm trying to test for brain-behavior correlations in some fMRI data. I have 34 clusters/ROIs that came from another analysis (independent, I hope, since it did not involve the behavioral measure to generate them).

I have 37 participants, 12 runs per participant (i.e., a run is a single session of neuroimaging data). For a few participants there are only 6 runs.

I have three experimental variables. Levels of variable A apply to an entire run, while levels of variables B and C apply to trials within a run. There are three levels of A and two of B and C.

I am trying to specify a model in which the contrast over levels in variable B in the brain data for a run is regressed on the same contrast in the behavioral data. I also want to know how this is modulated by levels of A and C.

I started with:

model <- lmer(B_contrast_brain ~ B_contrast_behavior*A*C + (1|subject/cluster), data=data)


I was getting some p-values that feel implausibly low to me, so I was trying to figure out what I'm doing wrong. One thing I realized: because I have two levels of C per run, run is also a grouping variable nested within subject. But I'm not sure of the correct way to specify that.

My best guess is:

  model <- lmer(B_contrast_brain ~ B_contrast_behavior*A*C + (1|subject/cluster) + (1|subject/run), data=data)


But I'm not sure. When I do that my p-values inch into the realm of plausibility, although a few of the effects are still very strong.

So my questions are:

1. Is the above the right specification?

2. Is there anything else that pops out that I'm doing wrong that could lead to spuriously strong effects?

Edit: I should add that since I asked this question I realized that I was getting the singular fit warning. Since paring down my fixed effects is not making that warning go away, this must be to do with how I specified my random effects. This might explain p-values that still feel implausibly low, but I guess I don't understand why I don't have enough data to have both cluster and run be random effects, when I have 34 clusters and 12 runs per subject, with cluster and run fully crossed.

Edit 2, per request:

Here is the output of summary for the first model. rel is what I have been calling B, type is A, and taught is C.

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: winsorized_rel_brain ~ rel_behavior * type * taught + (1 | subject/cluster)
Data: all_data_dx

REML criterion at convergence: 71101.8

Scaled residuals:
Min      1Q  Median      3Q     Max
-5.0352 -0.3851 -0.0614  0.3395  5.7847

Random effects:
Groups          Name        Variance Std.Dev.
cluster:subject (Intercept) 0.003915 0.06257
subject         (Intercept) 0.010999 0.10488
Residual                    0.771609 0.87841
Number of obs: 27466, groups:  cluster:subject, 1153; subject, 37

Fixed effects:
Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                 6.153e-02  2.170e-02  7.076e+01   2.836 0.005957 **
rel_behavior                4.071e-02  1.489e-02  2.706e+04   2.734 0.006264 **
typeN                      -8.669e-02  1.847e-02  2.712e+04  -4.693 2.70e-06 ***
typeY                      -1.113e-03  1.849e-02  2.707e+04  -0.060 0.951981
taughtU                    -1.387e-01  1.835e-02  2.706e+04  -7.556 4.29e-14 ***
rel_behavior:typeN         -5.066e-02  2.023e-02  2.708e+04  -2.504 0.012283 *
rel_behavior:typeY         -1.931e-02  2.231e-02  2.711e+04  -0.866 0.386571
rel_behavior:taughtU       -9.047e-02  2.336e-02  2.693e+04  -3.872 0.000108 ***
typeN:taughtU               2.093e-01  2.617e-02  2.707e+04   7.998 1.31e-15 ***
typeY:taughtU               7.923e-02  2.622e-02  2.707e+04   3.021 0.002518 **
rel_behavior:typeN:taughtU  7.945e-02  3.394e-02  2.695e+04   2.341 0.019227 *
rel_behavior:typeY:taughtU  1.676e-01  3.440e-02  2.711e+04   4.874 1.10e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) rl_bhv typeN  typeY  taghtU rl_b:N rl_b:Y rl_b:U typN:U typY:U r_:N:U
rel_behavir  0.028
typeN       -0.421 -0.034
typeY       -0.420 -0.034  0.493
taughtU     -0.423 -0.034  0.497  0.497
rl_bhvr:tyN -0.021 -0.734  0.101  0.023  0.025
rl_bhvr:tyY -0.020 -0.660  0.022  0.103  0.023  0.482
rl_bhvr:tgU -0.017 -0.649  0.021  0.023  0.018  0.468  0.434
typeN:tghtU  0.297  0.023 -0.705 -0.348 -0.701 -0.069 -0.016 -0.011
typeY:tghtU  0.296  0.025 -0.348 -0.705 -0.700 -0.018 -0.074 -0.014  0.491
rl_bhvr:N:U  0.011  0.449 -0.061 -0.014 -0.012 -0.607 -0.292 -0.683 -0.025  0.009
rl_bhvr:Y:U  0.013  0.428 -0.014 -0.067 -0.012 -0.310 -0.651 -0.670  0.008 -0.015  0.454


Here is the output for the second:

 Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: winsorized_rel_brain ~ rel_behavior * type * taught + (1 | subject/cluster) +      (1 | subject/run)
Data: all_data_dx

REML criterion at convergence: 68406.7

Scaled residuals:
Min      1Q  Median      3Q     Max
-7.4126 -0.3985 -0.0264  0.3608  6.9236

Random effects:
Groups          Name        Variance  Std.Dev.
cluster.subject (Intercept) 0.000e+00 0.0000000
run.subject     (Intercept) 1.091e-01 0.3303322
subject         (Intercept) 0.000e+00 0.0000000
subject.1       (Intercept) 1.638e-10 0.0000128
Residual                    6.784e-01 0.8236462
Number of obs: 27466, groups:  cluster:subject, 1153; run:subject, 443; subject, 37

Fixed effects:
Estimate Std. Error         df t value Pr(>|t|)
(Intercept)                 5.231e-02  2.976e-02  5.206e+02   1.758 0.079394 .
rel_behavior                4.339e-02  1.890e-02  1.869e+04   2.296 0.021689 *
typeN                      -9.501e-02  4.216e-02  5.238e+02  -2.253 0.024652 *
typeY                      -7.202e-03  4.223e-02  5.235e+02  -0.171 0.864643
taughtU                    -1.394e-01  1.722e-02  2.702e+04  -8.096 5.93e-16 ***
rel_behavior:typeN         -1.067e-01  2.563e-02  1.894e+04  -4.165 3.13e-05 ***
rel_behavior:typeY         -7.256e-02  2.801e-02  2.008e+04  -2.590 0.009597 **
rel_behavior:taughtU        4.737e-03  3.241e-02  1.452e+04   0.146 0.883789
typeN:taughtU               2.160e-01  2.465e-02  2.709e+04   8.761  < 2e-16 ***
typeY:taughtU               8.743e-02  2.475e-02  2.714e+04   3.532 0.000413 ***
rel_behavior:typeN:taughtU  5.267e-02  4.664e-02  1.504e+04   1.129 0.258884
rel_behavior:typeY:taughtU  1.129e-01  4.516e-02  1.769e+04   2.500 0.012416 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) rl_bhv typeN  typeY  taghtU rl_b:N rl_b:Y rl_b:U typN:U typY:U r_:N:U
rel_behavir  0.027
typeN       -0.706 -0.019
typeY       -0.705 -0.019  0.498
taughtU     -0.290 -0.045  0.205  0.204
rl_bhvr:tyN -0.020 -0.737  0.056  0.014  0.033
rl_bhvr:tyY -0.018 -0.675  0.013  0.056  0.030  0.498
rl_bhvr:tgU -0.020 -0.734  0.014  0.014  0.028  0.541  0.495
typeN:tghtU  0.202  0.031 -0.290 -0.143 -0.698 -0.076 -0.021 -0.019
typeY:tghtU  0.202  0.031 -0.142 -0.290 -0.695 -0.023 -0.095 -0.019  0.486
rl_bhvr:N:U  0.014  0.510 -0.039 -0.010 -0.019 -0.698 -0.344 -0.695 -0.025  0.013
rl_bhvr:Y:U  0.014  0.527 -0.010 -0.036 -0.020 -0.388 -0.691 -0.718  0.014 -0.016  0.499
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

• Please include the output of summary(model) for both models. Also please can you let us know how you determine whether a p-value is plausible ? Sep 17, 2021 at 20:41
• @RobertLong thanks for the engagement! I added both outputs. As for the p-values, if asked to weigh the likelihood that an interaction in my experimental variables (e.g. typeN:taughtU ) is determining a contrast in brain activation estimable with such precision, versus the likelihood that I'm making a mistake in reasoning somewhere, I think "I'm making a mistake" is more plausible. Sep 17, 2021 at 21:24
• Specifically, I think my mistake may be about the degrees of freedom in my denominator, but I can't figure out how to correctly specify them. A domain expert also pointed me to this blogpost, which may be affecting my inference as well: mumfordbrainstats.tumblr.com/post/126904300281/… Sep 17, 2021 at 21:27
• The 2nd model is singular. As for the denomnator degress of freedom, in a mixed model this is an approximation at best and woefully wrong at worst and underlines why I recommend not using p-values at all - they are little more than a test of sample size. Sep 18, 2021 at 8:16

From the description given it does sound as though you have repeated measures within run, which is not accounted for in the first model. However, the 2nd model:

B_contrast_brain ~ B_contrast_behavior*A*C + (1|subject/cluster) + (1|subject/run)


has the following features:

• fixed effects for B_contrast_behavior, A, C and all the interactions
• random intercepts for subject and also for cluster varying within levels of subject
• random intercepts for subject (again) and also for run varying within levels of subject.

So, you have specified random intercepts for subject twice, which may very well be the cause of the singular fit for the 2nd model. In order to just fit random intercepts for subject once you can use:

B_contrast_brain ~ B_contrast_behavior*A*C + (1|subject/cluster) + (1|subject:run)


where we specify random intercepts for run varying within levels of subject but without re-specifying random intercepts for subject again.

As for the p-values, in a mixed model these are approximations at best, so I would avoid them. Also, since you appear to have quite a large sample size, with rather small variance components for subject and cluster (based on the first model output) small p-values are to be expected.

• Thank you! Interestingly, when I run your suggested model with my winsorized DV (I didn't winsorize very aggressively; I brought the tails in to about 4 SD) I still got a singular fit that this time failed to converge, but when I ran it on my standardized DV with no winsorization it was fine. I was only winsorizing as an attempt to eliminate the possibility that (really large) deviations from normality of the residuals in the model were the source of my strong effects. Sep 18, 2021 at 12:27
• On reflection, the strongest effects here are expected, because they are not interactions with the behavior variable, but rather another way of representing contrasts that I also ran whole brain, and found effects in -- in different clusters, but still, if the three way interaction between A and B and C was powerful enough to generate a bunch of clusters that survived correction for multiple comparisons, it's not surprising that I see it reflected in a different way in my data here. Sep 18, 2021 at 12:32
• Do you have a link to where you describe your alternatives to p-values recommendations? Do you recommend confidence intervals? Some kind of bootstrappy simulation test? Sep 18, 2021 at 12:34
• I didn't recommend an alternative to p-values. I don't need to know the probability of observing my data if the null hypothesis is in fact true, and that's all a p-value tells me. Sep 18, 2021 at 13:50