Visualizing (generalized) linear mixed effects models, part 2 #rstats #lme4

In the first part on visualizing (generalized) linear mixed effects models, I showed examples of the new functions in the sjPlot package to visualize fixed and random effects (estimates and odds ratios) of (g)lmer results. Meanwhile, I added further features to the functions, which I like to introduce here. This posting is based on the online manual of the sjPlot package.

In this posting, I’d like to give examples for diagnostic and probability plots of odds ratios. The latter examples, of course, only refer to the sjp.glmer function (generalized mixed models). To reproduce these examples, you need the version 1.59 (or higher) of the package, which can be found at GitHub. A submission to CRAN is planned for the next days…

Fitting example models

The following examples are based on two fitted mixed models:

# fit model
library(lme4)
# create binary response
sleepstudy$Reaction.dicho <- sju.dicho(sleepstudy$Reaction, 
                                       dichBy = "md")
# fit first model
fit <- glmer(Reaction.dicho ~ Days + (Days | Subject),
             sleepstudy,
             family = binomial("logit"))

data(efc)
# create binary response
efc$hi_qol <- sju.dicho(efc$quol_5)
# prepare group variable
efc$grp = as.factor(efc$e15relat)
levels(x = efc$grp) <- sji.getValueLabels(efc$e15relat)
# data frame for 2nd fitted model
mydf <- na.omit(data.frame(hi_qol = as.factor(efc$hi_qol),
                           sex = as.factor(efc$c161sex),
                           c12hour = as.numeric(efc$c12hour),
                           neg_c_7 = as.numeric(efc$neg_c_7),
                           grp = efc$grp))
# fit 2nd model
fit2 <- glmer(hi_qol ~ sex + c12hour + neg_c_7 + (1|grp),
              data = mydf,
              family = binomial("logit"))

Summary fit1

Formula: Reaction.dicho ~ Days + (Days | Subject)
   Data: sleepstudy

     AIC      BIC   logLik deviance df.resid 
   158.7    174.7    -74.4    148.7      175 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.2406 -0.2726 -0.0198  0.2766  2.9705 

Random effects:
 Groups  Name        Variance Std.Dev. Corr 
 Subject (Intercept) 8.0287   2.8335        
         Days        0.2397   0.4896   -0.19
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.8159     1.1728  -3.254 0.001139 ** 
Days          0.8908     0.2347   3.796 0.000147 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
Days -0.694

Summary fit2

Formula: hi_qol ~ sex + c12hour + neg_c_7 + (1 | grp)
   Data: mydf

     AIC      BIC   logLik deviance df.resid 
  1065.3   1089.2   -527.6   1055.3      881 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7460 -0.8139 -0.2688  0.7706  6.6464 

Random effects:
 Groups Name        Variance Std.Dev.
 grp    (Intercept) 0.08676  0.2945  
Number of obs: 886, groups:  grp, 8

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.179036   0.333940   9.520  < 2e-16 ***
sex2        -0.545282   0.178974  -3.047  0.00231 ** 
c12hour     -0.005148   0.001720  -2.992  0.00277 ** 
neg_c_7     -0.219586   0.024108  -9.109  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
        (Intr) sex2   c12hor
sex2    -0.410              
c12hour -0.057 -0.048       
neg_c_7 -0.765 -0.009 -0.116

Diagnostic plots

Two new functions are added to both sjp.lmer and sjp.glmer, hence they apply to linear and generalized linear mixed models, fitted with the lme4 package. The examples only refer to the sjp.glmer function.

Currently, there are two type options to plot diagnostic plots: type = "fe.cor" to plot a correlation matrix between fixed effects and type = "re.qq" to plot a qq-plot of random effects.

Correlation matrix of fixed effects

To plot a correlation matrix of the fixed effects, use type = "fe.cor".

# plot fixed effects correlation matrix
sjp.glmer(fit2, type = "fe.cor")

qq-plot of random effects

Another diagnostic plot is the qq-plot for random effects. Use type = "re.qq" to plot random against standard quantiles. The dots should be plotted along the line.

# plot qq-plot of random effects
sjp.glmer(fit, type = "re.qq")

Probability curves of odds ratios

These plotting functions have been implemented to easier interprete odds ratios, especially for continuous covariates, by plotting the probabilities of predictors.

Probabilities of fixed effects

With type = "fe.pc" (or type = "fe.prob"), probability plots for each covariate can be plotted. These probabilties are based on the fixed effects intercept. One plot per covariate is plotted.

The model fit2 has one binary and two continuous covariates:

# plot probability curve of fixed effects
sjp.glmer(fit2, type = "fe.pc")

Probabilities of fixed effects depending on grouping level (random intercept)

With type = "ri.pc" (or type = "ri.prob"), probability plots for each covariate can be plotted, depending on the grouping level from the random intercept. Thus, for each covariate a plot for each grouping levels is plotted. Furthermore, with the show.se the standard error of probabilities can be shown. In this example, only the plot for one covariate is shown, not for all.

# plot probability curves for each covariate
# grouped by random intercepts
sjp.glmer(fit2,
          type = "ri.pc",
          show.se = TRUE)

Instead of faceting plots, all grouping levels can be shown in one plot:

# plot probability curves for each covariate
# grouped by random intercepts
sjp.glmer(fit2,
          type = "ri.pc",
          facet.grid = FALSE)

Outlook

These will be the new features for the next package update. For later updates, I’m also planning to plot interaction terms of (generalized) linear mixed models, similar to the existing function for visualizing interaction terms in linear models.