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"))
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
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
Two new functions are added to both
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
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.
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)
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)
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.