# Estimates radically change when including Random Slopes in Multiple Logistic Regression

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.

• Specifying random effects can be tricky, as different formulas allow for different potential correlations among estimates. Think carefully about exactly what you want to be modeling. See this answer on the lmer cheat sheet and pay attention to the difference between model M2 and model M3 in terms of how they handle correlations between intercept deviations and fixed-effect deviations for random effects. See this answer for enforcing independence of random-effect relations for 2 fixed effects. – EdM Feb 20 '19 at 17:27
• Also, from the number of observations (651) and the number of ID values (163), it seems that each ID only faced one set of the 4 combinations of Ambiguity and Uncertainty (with one missing data point). Your second model, with its large number of estimated random effects, might be overfitting. This difficulty here isn't with logistic regression per se, it's with what you want to and can adjust for with the mixed model. That choice depends on your understanding of the subject matter. – EdM Feb 20 '19 at 17:55
• For a start, I do not understand how can the intercept be ~10 when the data are only 0s and 1s. Am I missing something? – gsanroma Feb 20 '19 at 18:10
• @GerardSanroma That is a question about logistic regression per se: what you are modeling is the log of the odds ratio for PunishmentYes/No. That's what the logit link function is doing: it maps probabilities (from 0 to 1) onto the range of all real numbers. Thus intercepts and slopes can take on any real values. To convert from logit scale to probabilitiy, exponentiate to get the odds ratio, then convert the odds ratio to the corresponding probability. – EdM Feb 20 '19 at 19:57
• @EdM but still with a maximum slope of ~2 in the 1st model and the data in 0, 1, isn't the logit highly biased to one side of the probability (up to the effect of the random intercept) ? – gsanroma Feb 20 '19 at 20:15

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.

• EdM, see above some updated edits following your comments. You are completely right, introducing random slopes after including the interaction term cannot be performed due to the lack of a higher number of observations. So I will get satisfied by allowing my model to account for random intercepts! However, I still did not figure out what is going on with the estimation of my log odds... Any further suggestion on this? Thanks in advance! – Daniel Toribio Feb 25 '19 at 10:05
• @DanielToribio my sense is that you can describe your data quite well without including any random effects, just a standard logistic regression including the evidently important interaction term. For a start on keeping random effects if you still need to, see what happens if you simply omit the one case that did not complete all 4 combinations of Ambiguity and Uncertainty. The technical details behind optimization are beyond my expertise. See the comments above on your question from Ben Bolker, one of the lme4 authors, about other things to try and about posting your data. – EdM Feb 25 '19 at 14:21
• after agreeing with my coauthors, I shared a git repository with the relevant data. Feel free to have a look at it. – Daniel Toribio Feb 26 '19 at 9:18