Mixed Models 2

Setting Up

library(tidyverse)
library(lme4)
library(lmerTest)
library(afex)
library(MuMIn)

selfie <- read_csv('ruth_page_selfies.csv')

-- Column specification ----------------------------------------------------------------------------------------
cols(
  Angle = col_character(),
  UglyCat = col_character(),
  ID = col_character()
)

Logistic Mixed-Effects Regression

It’s generally necessary to convert text outcome variables to factors for generalized logistic regression.

selfie <- mutate(selfie, UglyCat = factor(UglyCat))

Let’s explore the data a little:

selfie %>% group_by(ID) %>% summarise(mean = mean(as.numeric(UglyCat)-1)) %>%
  ggplot(aes(x = mean))+
  geom_histogram(binwidth = .10)+
  theme_bw()

Clearly a lot of variation by person. Let’s see if we can figure out something about the differences in judging ‘ugly’ based on Angle across individuals…

selfie %>% group_by(ID, Angle) %>% 
  summarise(mean = mean(as.numeric(UglyCat)-1)) %>%
  pivot_wider(names_from = Angle, values_from = mean) %>%
  mutate(diff = FromBelow - Level) %>%
  ggplot(aes(x = diff))+
  geom_histogram(binwidth = .10)+
  theme_bw()

`summarise()` has grouped output by 'ID'. You can override using the `.groups` argument.

Also looks like quite a bit of variation in how Angle affects participants - for many there is little difference in the proportion of ‘ugly’ ratings assigned across the two levels of Angle, but for many participants more ratings of ‘ugly’ were given when the Angle was ‘FromBelow’.

Now to run a model:

ugly_mdl <- glmer(UglyCat ~ Angle + (1+Angle|ID),
                  data = selfie,
                  family = 'binomial')

Model failed to converge with max|grad| = 0.175575 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

This model doesn’t converge. Let’s take a peek at the summary…

summary(ugly_mdl)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie

     AIC      BIC   logLik deviance df.resid 
  1712.4   1739.2   -851.2   1702.4     1562 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8878 -0.5480  0.1709  0.5498  3.1089 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 ID     (Intercept) 3.683    1.919         
        AngleLevel  1.922    1.386    -0.39
Number of obs: 1567, groups:  ID, 98

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.9195887  0.0006310    1457   <2e-16 ***
AngleLevel  -1.5042298  0.0006313   -2383   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr)
AngleLevel -0.001
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.175575 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

It’s not immediately clear what’s wrong, but we should be very cautious about non-convergence.

We’ll try some optimizers to see if we can get a fit using different estimation strategies.

all_fit(ugly_mdl)

bobyqa. : [OK]
Nelder_Mead. : 

Model failed to converge with max|grad| = 0.112066 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

[OK]
optimx.nlminb : 

unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues

[OK]
optimx.L-BFGS-B : [OK]
nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
nloptwrap.NLOPT_LN_BOBYQA : [OK]
nmkbw. : [ERROR]
$bobyqa.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4533 1739.2379 -851.2267 1702.4533      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.900         
        AngleLevel  1.375    -0.38
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9049      -1.4834  

$Nelder_Mead.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4213 1739.2059 -851.2106 1702.4213      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.919         
        AngleLevel  1.386    -0.39
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9194      -1.5042  
optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 

$optimx.nlminb
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4250 1739.2096 -851.2125 1702.4250      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.929         
        AngleLevel  1.387    -0.40
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9164      -1.4946  
optimizer (optimx) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 

$`optimx.L-BFGS-B`
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4533 1739.2379 -851.2267 1702.4533      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.900         
        AngleLevel  1.375    -0.38
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9049      -1.4834  

$nloptwrap.NLOPT_LN_NELDERMEAD
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4533 1739.2379 -851.2267 1702.4533      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.900         
        AngleLevel  1.375    -0.38
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9049      -1.4834  

$nloptwrap.NLOPT_LN_BOBYQA
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
      AIC       BIC    logLik  deviance  df.resid 
1712.4533 1739.2379 -851.2267 1702.4533      1562 
Random effects:
 Groups Name        Std.Dev. Corr 
 ID     (Intercept) 1.900         
        AngleLevel  1.375    -0.38
Number of obs: 1567, groups:  ID, 98
Fixed Effects:
(Intercept)   AngleLevel  
     0.9049      -1.4834  

$nmkbw.
<packageNotFoundError in loadNamespace(name): there is no package called ‘dfoptim’>

This wall of output shows that (a) several optimizers can run the model and (b) nothing really changes too much when using different optimizers (i.e., fixed and random effect coefficients are very similar). This is good news.

Re-run the model with the bobyqa optimizer.

ugly_mdl <- glmer(UglyCat ~ Angle + (1+Angle|ID),
                  data = selfie,
                  family = 'binomial',
                  control = glmerControl(optimizer = 'bobyqa'))

summary(ugly_mdl)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: UglyCat ~ Angle + (1 + Angle | ID)
   Data: selfie
Control: glmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  1712.5   1739.2   -851.2   1702.5     1562 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8757 -0.5489  0.1733  0.5502  3.0998 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 ID     (Intercept) 3.61     1.900         
        AngleLevel  1.89     1.375    -0.38
Number of obs: 1567, groups:  ID, 98

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.9049     0.2199   4.115 3.87e-05 ***
AngleLevel   -1.4834     0.2001  -7.414 1.22e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr)
AngleLevel -0.486

Let’s take a peek at the variance explained by the model:

r.squaredGLMM(ugly_mdl)

The null model is correct only if all variables used by the original model remain unchanged.

                   R2m       R2c
theoretical 0.07445407 0.5550465
delta       0.06778359 0.5053189

The ‘delta’ value is probably the better one to report, as it is based more closely on the observed data.

Let’s work on effect sizes:

fixefs <- cbind(fixef(ugly_mdl), confint(ugly_mdl, parm = "beta_"))

Computing profile confidence intervals ...

fixefs

                            2.5 %    97.5 %
(Intercept)  0.9049237  0.4773409  1.380091
AngleLevel  -1.4834319 -1.9211854 -1.093406

Now exponentiate to get odds ratios…

exp(fixefs)

                          2.5 %    97.5 %
(Intercept) 2.4717433 1.6117827 3.9752646
AngleLevel  0.2268578 0.1464333 0.3350733