How to decide whether to include a random interaction for the random effect in a Linear Mixed Model?

Question

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

Answer 1

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.