1
$\begingroup$

I am trying to understand whether an experimental variable that we will call "candy_position" has an effect on the reaction times (called here "rt") of people during a task. All of the participants have been exposed to the 2 levels of the candy_position variable multiple times, every time recording the rt.

rt = reaction times, continuos variable; candy_position = categorical variable with two levels; subject = categorical variable representing 50 different people that participated in the study.

I am currently analyzing the results in R with lme4. I could specify the model in two ways:

lmer(rt ~ candy_position + (1 | subject), data = my_data,  REML = F) # model A
lmer(rt ~ candy_position + (1 | subject)+ (1 | subject:candy_position), data = my_data,  REML = F) # model B

"model A" has a random intercept for every subject, accounting for the fact that every subject could have a different mean rt and, in a way, preventing the effect of candy_position on rt to be masked by the variability of reaction times between subjects.

"Model B" is "model A" plus a random intercept for every subject interacting with candy_position. This means that we are also capturing the different ways in which candy_position has an effect on every single subject's reaction times.

"Model A" and "model B" have different results for my fixed effect (candy_position associated p-value is <0.05 with "model A" but >0.05 with "model B"). Is it a matter of degrees of freedom being more for the fixed effect in "model A" than "model B"?

More importantly, I don't know which one I should use. On one hand, I am inclined to use "model A", because I am interested in finding whether on average candy_position has an effect on the reaction times of people, not whether this is true more for some people and less, or not at all true, for others. On the other hand, I don't see how adding the random interaction (1 | subject:candy_position) should diminish the amount of variance explained by the candy_position fixed term.

Thank you for any possible explanation/advice you could give me. (clearly) I am not a statistician, I will be even more grateful if you will consider that while replying.

Model A summary:

Linear mixed model fit by maximum likelihood. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: rt ~ candy_position + (1 | subject)
   Data: my_data

      AIC       BIC    logLik  deviance  df.resid 
 205401.2  205432.0 -102696.6  205393.2     16145 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8898 -0.6375 -0.1568  0.4822  6.2386 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 46959    216.7   
 Residual             19292    138.9   
Number of obs: 16149, groups:  subject, 30

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           673.663     39.587    30.023   17.02  < 2e-16 ***
candy_positionvalid    -6.826      2.314 16119.060   -2.95  0.00318 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
cndy_pstnvl -0.020

Model B summary:

Linear mixed model fit by maximum likelihood. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: rt ~ candy_position + (1 | subject) + (1 | subject:candy_position)
   Data: my_data

      AIC       BIC    logLik  deviance  df.resid 
 205174.4  205212.9 -102582.2  205164.4     16144 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9957 -0.6334 -0.1572  0.4751  6.0579 

Random effects:
 Groups                 Name        Variance Std.Dev.
 subject:candy_position (Intercept)   799.9   28.28  
 subject                (Intercept) 47148.6  217.14  
 Residual                           18934.3  137.60  
Number of obs: 16149, groups:  subject:candy_position, 60; subject, 30

Fixed effects:
                    Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)          673.712     40.001  30.519  16.842   <2e-16 ***
candy_positionvalid   -5.403      7.657  30.007  -0.706    0.486    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
cndy_pstnvl -0.093
$\endgroup$
1
  • $\begingroup$ This is not a programming question, and is better suited for statsexchange (crossvalidation). Saying that testing random effects is in general difficult. You can't base your decision on the fixed effect, but instead should test the random effect directly. Disregarding bayesian methods, several approximations do exist but the number of implementations are sparse. The simplest (although time-consuming) method is to just use bootstrapping methods. Several examples exist in bootMer function for how to perform such a test. $\endgroup$
    – Oliver
    May 31, 2021 at 17:01

1 Answer 1

1
$\begingroup$

The first model makes sense.

The second model has a number of issues. First, it is highly questionable to include a fixed effect, candy_position in this case - which is your main exposure for which you are interested in statistical inference, in an interaction which is used as a grouping variable for random intercepts. This will result in variance that was otherwise attributed to the fixed part of the model, also being shared by the random part, which is rarely what you want. Second, To specify (1 | subject) + (1 | subject:candy_position) you are saying that each level of candy_position of which there seems to be only 2, "belongs" to one and only one level of subject - that is, that candy_position is nested within subject. Note that (1 | subject) + (1 | subject:candy_position) is exactly the same as (1 | subject/candy_position). So that does not seem to be the case here.

$\endgroup$
2
  • $\begingroup$ Thank you very much for your answer @RobertLong, it was really helpful! I am puzzled by the fact that the results of model A ( rt ~ candy_position + (1 | subject) ) become similar to the corresponding paired t.test only if I use the "aggregated" data (mean by subject and condition). When I run the model on the "disaggregated" data ( multiple data points for every subject for every condition) the results significantly differ (even the t values). Do you think it is better either to edit the question to add this or to open a new one? $\endgroup$ Jun 2, 2021 at 10:29
  • $\begingroup$ You're welcome :) I would definitely ask a new question about that. $\endgroup$ Jun 2, 2021 at 10:58

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.