Generalized Linear Mixed Model: Random slope inverts effect

Question

I struggle with the analysis of my very skewed data with linear mixed models in R. Since the original data is for actual research, I can't share it with you, but I have created a fake dataset, that resembles the distribution of my data: Let's assume, we give 1000 amateur dart players 4 throws and measure, if they can hit the board. Then, they have to drink 3 shots and try again (so we have 8 throws, 4 sober, 4 drunk). This would result in a dataset like this:

EDIT: Forgot to include the libraries

library(tidyverse)
library(lmerTest)

set.seed(1)

size <- 1000

dart_wide <- tibble(
  ID = factor(paste0('ID_', seq(size))),
  sober = sample(c(.75, 1), size = size, replace = TRUE, prob = c(.1, .9)),
  drunk = sapply(sober, simplify = TRUE, FUN = function(x) {
    if (x == 0.75) sample(c(0, .25, .75, 1), size = 1, prob = c(.01, .01, .02, .06))
    else if (x == 1) sample(c(0, .25, .5, .75, 1), size = 1, prob = c(.05, .05, .05, .15, .6))
  }),
  difference = drunk - sober
) 
  
dart <- pivot_longer(dart_wide, sober:drunk, names_to = 'condition', values_to = 'score')
dart <- dart %>% rbind(dart) %>% rbind(dart) %>% rbind(dart) %>%
  group_by(ID, condition) %>%
  mutate(
    condition = factor(condition, levels = c('sober', 'drunk')),
    throw = seq(n()),
    hit = as.numeric(score >= (throw/4))
  ) %>%
  ungroup %>%
  select(ID, condition, throw, hit) %>%
  arrange(ID, condition, throw)
dart

> # A tibble: 8,000 x 4
>    ID    condition throw   hit
>    <fct> <fct>     <int> <dbl>
>  1 ID_1  sober         1     1
>  2 ID_1  sober         2     1
>  3 ID_1  sober         3     1
>  4 ID_1  sober         4     1
>  5 ID_1  drunk         1     1
>  6 ID_1  drunk         2     1
>  7 ID_1  drunk         3     1
>  8 ID_1  drunk         4     1
>  9 ID_10 sober         1     1
> 10 ID_10 sober         2     1
> # ... with 7,990 more rows

In which most of the players perform worse in the "drunk" condition, compared to the "sober" condition

mean(dart_wide$difference)

> [1] -0.146

However the distribution does not even resemble a normal distribution:

dart %>% 
  group_by(ID, condition) %>% 
  summarize(accuracy = mean(hit), .groups = 'drop') %>% 
  ggplot(aes(accuracy)) +
  facet_wrap(~condition) +
  geom_bar()

If I now analyze the data with a generalized linear mixed model using only a random intercept, I get the expected results: The accuracy of the participants was lower when drunk (β = -2.4982):

dart_model_intercept <- glmer(
  hit ~ condition + (1|ID), 
  data = dart,
  family = binomial(link = 'logit')
)
summary(dart_model_intercept)



> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
>  Family: binomial  ( logit )
> Formula: hit ~ condition + (1 | ID)
>    Data: dart
> 
>      AIC      BIC   logLik deviance df.resid 
>   4217.8   4238.7  -2105.9   4211.8     7997 
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -7.7350  0.0632  0.1999  0.2203  1.3237 
> 
> Random effects:
>  Groups Name        Variance Std.Dev.
>  ID     (Intercept) 3.117    1.765   
> Number of obs: 8000, groups:  ID, 1000
> 
> Fixed effects:
>                Estimate Std. Error z value Pr(>|z|)    
> (Intercept)      4.8973     0.1597   30.67   <2e-16 ***
> conditiondrunk  -2.4982     0.1221  -20.46   <2e-16 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> Correlation of Fixed Effects:
>             (Intr)
> conditndrnk -0.763

Okay, great. HOWEVER, if I include the random slope for condition, the model estimates a contra-intuitive result: apparently, participants perform better when drunk than when sober (β = 0.6111)...

dart_model_slope <- glmer(
  hit ~ condition + (condition|ID), 
  data = dart,
  family = binomial(link = 'logit')
)
summary(dart_model_slope)

> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
>  Family: binomial  ( logit )
> Formula: hit ~ condition + (condition | ID)
>    Data: dart
> 
>      AIC      BIC   logLik deviance df.resid 
>   3963.2   3998.1  -1976.6   3953.2     7995 
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -5.9729  0.0952  0.1672  0.1674  1.4121 
> 
> Random effects:
>  Groups Name           Variance  Std.Dev. Corr
>  ID     (Intercept)     0.004656 0.06824      
>         conditiondrunk 13.328751 3.65086  1.00
> Number of obs: 8000, groups:  ID, 1000
> 
> Fixed effects:
>                Estimate Std. Error z value Pr(>|z|)    
> (Intercept)      3.5679     0.0970  36.782   <2e-16 ***
> conditiondrunk   0.6111     0.4150   1.472    0.141    
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> Correlation of Fixed Effects:
>             (Intr)
> conditndrnk -0.234

However, this is not the case, as we have seen that overall, participants perform worse (direct comparison of performance using mean())

Apparently, the model estimates positive slopes for participants, that don't change their performance between conditions

dart_wide_with_slope <- coef(dart_model_slope)$ID %>% 
  rownames_to_column('ID') %>% 
  select(ID, random_slope = conditiondrunk) %>% 
  left_join(dart_wide %>% select(ID, difference), by = 'ID')
dart_wide_with_slope %>% as_tibble

> # A tibble: 1,000 x 3
>    ID      random_slope difference
>    <chr>          <dbl>      <dbl>
>  1 ID_1           1.13        0   
>  2 ID_10          1.13        0   
>  3 ID_100        -4.15       -0.75
>  4 ID_1000        1.13        0   
>  5 ID_101        -2.12       -0.25
>  6 ID_102         1.13        0   
>  7 ID_103         1.13        0   
>  8 ID_104         0.964       0.25
>  9 ID_105         1.13        0   
> 10 ID_106         1.13        0   
> # ... with 990 more rows

So this then results in a positive overall weight for condition: drunk. But this is just objectively not the case/correct.

Now finally I reached my question:

1. From my understanding, random slopes are calculated by distributing the performance around the mean (super simplified). Is there a way to prevent this? This should probably resolve this issue
2. Is there a better analysis I can use (maybe use a different link in the binomial function?
3. Does it even make sense to use a GLMM for this distribution of data and if not, what would you propose?

Answer 1

Interesting problem! It seems to me that you would be better off to treat your response as a binomial rather than a binary response.

A binomial response would imply that, for each subject, you keep track of the number of hits (out of the number of draws) for each condition. For example:

ID    condition    hits    draws
1       sober       3        4
1       drunk       1        4 
2       sober       4        4
2       drunk       0        4
etc.

The dataset for the binomial response will have only two rows for each subject: one for the sober condition and one for the drunk condition. The number of draws for each condition (i.e., 4) will be recorded under the draws column. The number of hits (out of 4 draws) for each condition will be recorded under the hits column.

To fit a mixed effects binomial logit model, you would need a syntax of the form:

mod1 <- glmer(cbind(hits, draws - hits) ~ 1 + condition + (1|ID), 
            family = binomial(link="logit"), 
            data = dart) 

or

mod2 <- glmer(cbind(hits, draws - hits) ~ 1 + condition + (1 +  condition|ID), 
            family = binomial(link="logit"), 
            data = dart) 

Don’t forget you are dealing with a model which provides a subject-specific interpretation for your estimated effect. In other words, you are looking at odds of getting a hit in any of 4 trials for the typical subject (i.e., the subject with a random intercept of 0) and comparing those odds between the drunk and sover conditions for that subject. These are so-called conditional odds. To get marginal odds, you can use the mixed_model() function from the GLMMadaptive package of R. Something like this:

 Mod1 <- mixed_model(fixed = cbind(hits, draws - hits) ~ condition, 
                     random = ~ 1 | ID,  
                     data = darts,
                     family = binomial())

 
 Mod2 <- mixed_model(fixed = cbind(hits, draws - hits) ~ condition, 
                     random = ~ condition | ID,  
                     data = darts,
                     family = binomial())

Then you can compare the two models with the anova() function via a likelihood ratio test to see which one is more appropriate for your data:

anova(Mod1, Mod2)

According to https://drizopoulos.github.io/GLMMadaptive/articles/GLMMadaptive_basics.html, "the likelihood ratio test between the two models is computed with function anova(). When two "MixMod" objects are provided, the function assumes that the first object represents the model under the null hypothesis, and the second object the model under the alternative". Thus, in your case, a "small" p-value (e.g., smaller than 0.05) would favor the more complex Mod2. Note that, if you wanted no correlation between the random intercepts and random slopes of condition in your model, you would specify that model by using the || symbol in the specification of the random argument:

Mod3 <- mixed_model(fixed = cbind(hits, draws - hits) ~ condition, 
                     random = ~ condition || ID,  
                     data = darts,
                     family = binomial())

See https://drizopoulos.github.io/GLMMadaptive/ - what you want for marginal odds ratios is the marginal_coefs() function in the GLMMadaptive package.

Say your data indeed support Mod2 - then you would apply the marginal_coefs() function like this:

marginal_coefs(Mod2, std_errors = TRUE)

These marginal coefficients are reported on the log odds scale - to get marginal odds ratios, you will need to exponentiate them.

The bottom line is that you shouldn’t panic if the results you get from your model for the typical subject look a certain way, especially if your interest is in seeing what happens across subjects rather than for a typical subject.

So you will need to figure out whether you are interested in reporting a marginal odds ratio and/or a conditional odds ratio contrasting the odds of a hit in 4 draws among your drunk and sober conditions. Also, you will need to question your own assumptions: Are your subjects champion dart throwers who would do really well when sober but terribly when drunk? Or are they people who have never thrown darts before in their life and might do better when drunk than when sober? Your model does not really capture their prior level of know how (e.g., years of experience throwing darts).

There may be a learning effect at play here as well - subjects may get better at throwing darts just because they keep throwing. Would they still get better after drinking (but not as much as they would without drinking)? It would have been interesting to investigate that by having some study subjects not drink before their second set of 4 throws.