2
$\begingroup$

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
Improve this question
$\endgroup$
4
  • $\begingroup$ 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 ? $\endgroup$
    – Robert Long
    Commented Sep 17, 2021 at 20:41
  • $\begingroup$ @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. $\endgroup$
    – Katie
    Commented Sep 17, 2021 at 21:24
  • $\begingroup$ 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/… $\endgroup$
    – Katie
    Commented Sep 17, 2021 at 21:27
  • $\begingroup$ 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. $\endgroup$
    – Robert Long
    Commented Sep 18, 2021 at 8:16

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Katie
    Commented Sep 18, 2021 at 12:27
  • $\begingroup$ 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. $\endgroup$
    – Katie
    Commented Sep 18, 2021 at 12:32
  • $\begingroup$ 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? $\endgroup$
    – Katie
    Commented Sep 18, 2021 at 12:34
  • $\begingroup$ 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. $\endgroup$
    – Robert Long
    Commented Sep 18, 2021 at 13:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.