Estimates radically change when including Random Slopes in Multiple Logistic Regression

Question

I am examining the fixed effects of two within-subject experimental manipulations (i.e., Ambiguity 0 = No / 1 = Yes, and Uncertainty 0 = No / 1 = Yes) on a dichotomized variable (i.e., Punishment, 0 = No / 1 = Yes) through a multiple logistic regression. Since the design includes repeated measures, I consider my subject's ID as grouping factor.

My variable Punishment was dichotomized as follows:

##Example for one of the repeated measures:
data1$Punishment.R1 <- ifelse(data1$R1>0,1,ifelse(is.na(data1$R1),NA,0))

##Subset of resultant dataframe:
library(dplyr)
subset_data1 <- data1[order(data1$ID),] %>% select(c("ID","R1","Punishment.R1"))%>%head(20)
subset_data1

##    ID R1 Punishment.R1
##145  1 15             1
##146  2 14             1
##147  3  0             0
##148  4 15             1
##149  5 10             1
##150  6 15             1
##151  7 18             1
##152  8 12             1
##153  9  1             1
##154 10 15             1
##155 11  6             1
##156 12 13             1
##157 13 15             1
##158 14 12             1
##159 15  0             0
##160 16  0             0
##161 17 15             1
##162 18 15             1
##163 19  3             1
##164 20 15             1

These are the frequencies tables for the 4 different within-subject conditions:

#NO Ambiguity / NO Uncertainty
table(data1$Punishment.R1)/nrow(data1)
##         0         1 
## 0.1463415 0.8475610

#NO Ambiguity / YES Uncertainty
table(data1$Punishment.R2)/nrow(data1)
##         0         1 
## 0.1585366 0.8292683

#YES Ambiguity / NO Uncertainty 
table(data1$Punishment.R3)/nrow(data1)
##         0         1 
## 0.2682927 0.7256098

#YES Ambiguity / YES Uncertainty
table(data1$Punishment.R4)/nrow(data1)
##         0         1 
## 0.2195122 0.7743902

After reshaping my data from wide to long format and creating my two within-subject factors, this is how my dataframe dichm used in following analyses looks like for the first 5 subjects:

library(dplyr)
subset_dichm <- dichm[order(dichm$ID),] %>% select(c("ID","Round","Ambiguity","Uncertainty","Punishment")) %>% head(20)
subset_dichm

    ID         Round Ambiguity Uncertainty Punishment
52   1 Punishment.R1         0           0          1
234  1 Punishment.R2         0           1          1
416  1 Punishment.R3         1           0          1
598  1 Punishment.R4         1           1          1
53   2 Punishment.R1         0           0          1
236  2 Punishment.R2         0           1          1
417  2 Punishment.R3         1           0          1
599  2 Punishment.R4         1           1          1
54   3 Punishment.R1         0           0          0
237  3 Punishment.R2         0           1          0
418  3 Punishment.R3         1           0          0
600  3 Punishment.R4         1           1          0
55   4 Punishment.R1         0           0          1
238  4 Punishment.R2         0           1          1
419  4 Punishment.R3         1           0          1
602  4 Punishment.R4         1           1          0
56   5 Punishment.R1         0           0          1
239  5 Punishment.R2         0           1          1
420  5 Punishment.R3         1           0          1
603  5 Punishment.R4         1           1          1

All the analyses presented below were run in R, using the lme4 package (version 1.1.17).

My first step was to allow random intercepts into my model:

H1.RI <- glmer(Punishment~Ambiguity + Uncertainty + (1|ID), data=dichm,family = binomial(link = logit))

##Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##Family: binomial  ( logit )
##Formula: Punishment ~ Ambiguity + Uncertainty + (1 | ID)
##Data: dichm
##AIC      BIC   logLik deviance df.resid 
##395.5    413.4   -193.7    387.5      647 
##Scaled residuals: 
##Min      1Q  Median      3Q     Max 
##-4.4905  0.0051  0.0063  0.0213  1.9409 
##Random effects:
## Groups Name        Variance Std.Dev.
## ID     (Intercept) 125      11.18   
##Number of obs: 651, groups:  ID, 163
##Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
##(Intercept)   10.0213     1.0292   9.737  < 2e-16 ***
##Ambiguity1    -2.4270     0.4917  -4.936 7.98e-07 ***
##Uncertainty1   0.4281     0.4131   1.036      0.3    
##---
##Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##Correlation of Fixed Effects:
##            (Intr) Ambgt1
##Ambiguity1  -0.602       
##Uncertanty1 -0.084 -0.073

As you can see, I observed a significant negative effect of the Ambiguity manipulation on Punishment.

The second model I tested included the effects of both experimental manipulations as random effects to account for potential intra-individual response patterns. See output below:

H1.RS <- glmer(Punishment~Ambiguity + Uncertainty + (1+Uncertainty+Ambiguity|ID), data=dichm,family = binomial(link = logit))

##Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##Family: binomial  ( logit )
##Formula: Punishment ~ Ambiguity + Uncertainty + (1 + Uncertainty + Ambiguity |      ID)
##   Data: dichm
## AIC      BIC   logLik deviance df.resid 
## 324.0    364.3   -153.0    306.0      642 
##Scaled residuals: 
##Min      1Q  Median      3Q     Max 
##-4.4929  0.0000  0.0000  0.0096  0.7864 
##Random effects:
## Groups Name         Variance Std.Dev. Corr     
## ID     (Intercept)    92.14   9.599            
##        Uncertainty1  964.71  31.060   1.00     
##        Ambiguity1   4833.49  69.523   0.36 0.36
##Number of obs: 651, groups:  ID, 163
##Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
##(Intercept)     8.042      1.278   6.295 3.08e-10 ***
##Ambiguity1      4.912      1.739   2.825  0.00473 ** 
##Uncertainty1   18.325      2.570   7.131 9.96e-13 ***
##---
##Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##Correlation of Fixed Effects:
##            (Intr) Ambgt1
##Ambiguity1  -0.637       
##Uncertanty1  0.252 -0.110

Surprisingly, there is not only a substantial change in the size of my estimates, but a complete reversal of the directionality of the Ambiguity effect (from -2.4270 to 4.912).

I am struggling to understand why is this the case, and how should I proceed from here in order to clarify what is happening with my data. I conducted a model comparison which seems to point at the random slope model as the better fit to the data:

anova(H1.RI,H1.RS)

##H1.RI: Punishment ~ Ambiguity + Uncertainty + (1 | ID)
##H1.RS: Punishment ~ Ambiguity + Uncertainty + (1 + Uncertainty + Ambiguity | ID)
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
##H1.RI   4 395.45 413.37 -193.73   387.45                            
##H1.RS  9 323.99 364.30 -153.00   305.99 81.462      5  4.149e-16 ***
##---
##Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I am not really familiar with logistic regression so I might be missing something specific from this kind of analysis that applies to the test of mixed models. Also, I have read in other threads that it could be a case of a Simpson's Paradox (see here), but I would like to know if there is any other plausible explanation given the data.

Any comment, reference, or redirection to a similar thread is highly appreciated.

Data can be accessed here: https://github.com/toribio-florez/Dichm-Data-Testing

---

To the recommendations added by @EdM:

Indeed, the glm gives your expected log odds estimates, so we can be more certain about code glitches not being the problem here:

H1.glm <- glm(Punishment~Ambiguity + Uncertainty, data=dichm, family = binomial(link = logit))
summary(H1.glm)

##Deviance Residuals: 
##    Min       1Q   Median       3Q      Max  
##-1.9603   0.5627   0.5930   0.7317   0.7690  
##Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
##(Intercept)    1.6491     0.1805   9.135  < 2e-16 ***
##Ambiguity1    -0.5819     0.2006  -2.901  0.00371 ** 
##Uncertainty1   0.1139     0.1976   0.577  0.56408    
##---
##Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##(Dispersion parameter for binomial family taken to be 1)
##Null deviance: 650.97  on 650  degrees of freedom
##Residual deviance: 642.04  on 648  degrees of freedom
##(5 observations deleted due to missingness)
##AIC: 648.04
##Number of Fisher Scoring iterations: 4

I also followed your recommendations regarding using different optimizers used by, in this case, glmer, but they do not seem to substantially change the estimates of the model:

H1.RI <- glmer(Punishment~Ambiguity + Uncertainty + (1|ID), data=dichm,family = binomial(link = logit))
test1 <- update(H1.RI, control=glmerControl(optimizer="nloptwrap"))

##             Estimate Std. Error z value Pr(>|z|)    
##(Intercept)   10.0215     1.0287   9.742  < 2e-16 ***
##Ambiguity1    -2.4271     0.4916  -4.937 7.95e-07 ***
##Uncertainty1   0.4281     0.4131   1.036      0.3    

test2 <- update(H1.RI, control=glmerControl(optimizer="bobyqa"))

##             Estimate Std. Error z value Pr(>|z|)    
##(Intercept)   10.0213     1.0278   9.751  < 2e-16 ***
##Ambiguity1    -2.4270     0.4915  -4.937 7.92e-07 ***
##Uncertainty1   0.4281     0.4131   1.036      0.3    

Perhaps for glmer it is necessary to specify two optimizers, let me know if this is the case.

Finally, we specified our model following our pre-registered hypotheses regarding the independent effects of our manipulations, but it is true that this does not seem to be the case in our data, since the Ambiguity:Uncertainty interaction term reaches statistical significance:

H1.RI <- glmer(Punishment~Ambiguity + Uncertainty + Ambiguity*Uncertainty + (1|ID), data=dichm,family = binomial(link = logit))

                        Estimate Std. Error z value Pr(>|z|)    
##(Intercept)              10.8506     1.1500   9.436  < 2e-16 ***
##Ambiguity1               -3.3916     0.7153  -4.742 2.12e-06 ***
##Uncertainty1             -0.6084     0.6516  -0.934   0.3505
##Ambiguity1:Uncertainty1   1.8077     0.8808   2.052   0.0401 * 

But when introducing random slopes, I receive the expected error: Error: number of observations (=651) < number of random effects (=652) for term (1 + Uncertainty + Ambiguity + Ambiguity * Uncertainty | ID); the random-effects parameters are probably unidentifiable.

Answer 1

There are three issues here, one technical and two statistical.

The technical issue is that lmer seems to be giving incorrect values for the fixed effects. Based on your frequency tables, the log-odds in each combination of Ambiguity and Uncertainty is:

              Ambiguity
              No    Yes
Uncertainty
No            1.75  1.00
Yes           1.65  1.26

Logistic regression models the log-odds, so the fixed-effects intercept (estimate for no Ambiguity, no Uncertainty) should be close to 1.75. Your mixed models are providing extremely and unrealistically high values of that intercept (10 or 8 on the log-odds scale, with corresponding odds in the thousands), as @GerardSanroma noted. At first I wondered whether your use of numeric ID codes for the individuals might be involved, but it seems that lmer() converts those to a factor. You might nevertheless want to double-check by explicitly converting the ID values to a factor. It's usually safest to specify yourself the numeric/factor/ordinal/logical classes of predictors, rather than count on software to make the conversion for you.

To rule out hidden glitches in your data or code (I didn't find any), fit a standard fixed-effect glm() model to your data (ignoring the random effects for now). You should get a value about 1.7 for the intercept, a small coefficient for Uncertainty (about 0.1) and a larger coefficient for Ambiguity (about 0.6).

I think, however, we have to look elsewhere for the technical problem with your mixed models. The default optimization in lmer() can find a local rather than a global optimum. This phenomenon is illustrated on this page and this page. As noted on those pages, if that's the case you can try specifying a different optimizer via a parameter setting within the call to lmer():

control=lmerControl(optimizer="nloptwrap")

The first statistical issue is one that you didn't raise but that is seen in the table of log-odds above. Your model for the fixed effects included only direct effects for Ambiguity and Uncertainty. The implicit assumption is that the effect of Ambiguity is independent of the level of Uncertainty, and vice-versa. That doesn't seem to be the case in your results. On the log-odds scale (which the logistic model is using), Ambiguity has almost twice the magnitude of effect in Uncertainty=No as in Uncertainty=Yes. The direction of the effect of Uncertainty differs depending on Ambiguity. It seems important to allow for and to test an interaction term (Ambiguity:Uncertainty), even more important than to allow for random effects. Leaving out an important interaction term can easily lead to Simpson's paradox.

The second statistical issue is what you can reasonably hope to control for with random effects. With an Ambiguity:Uncertainty fixed-effect interaction term that seems likely to be significant, it would be unwise to try to allow for random main slopes without also including corresponding random interactions. Yet you have at most 1 of each of the 4 Ambiguity-Uncertainty combinations tested per individual. You have fewer observations than the number of effects that you are trying to estimate: the fixed effects of Ambiguity, Uncertainty and their interaction, and corresponding random effects for each individual. Thus there is no hope of "accounting for the intra-individual sensitivity to these manipulations." With your first model (intercept only) you could allow for differences among individuals in their overall tendency to choose Punishment=yes, which might provide some advantage. It doesn't seem wise to go much beyond that in terms of random effects with your data.