Introduction: Marginal Effects for Random Effects Models

Daniel Lüdecke

2019-07-01

This vignette shows how to calculate marginal effects that take the random-effect variances for mixed models into account.

Marginal effects for mixed effects models

Basically, the type of predictions, i.e. whether to account for the uncertainty of random effects or not, can be set with the type-argument. The default, type = "fe", means that predictions are on the population-level and do not account for the random effect variances. Intervals are confidence intervals for the predicted values.

library(ggeffects)
library(lme4)
data(sleepstudy)
m <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)

pr <- ggpredict(m, "Days")
pr
#> 
#> # Predicted values of Reaction
#> # x = Days
#> 
#>  x predicted std.error conf.low conf.high
#>  0   251.405     6.824  238.031   264.779
#>  1   261.872     6.786  248.572   275.173
#>  2   272.340     7.094  258.436   286.243
#>  3   282.807     7.705  267.705   297.909
#>  5   303.742     9.582  284.962   322.521
#>  6   314.209    10.733  293.173   335.245
#>  7   324.676    11.974  301.208   348.144
#>  9   345.611    14.630  316.936   374.285
#> 
#> Adjusted for:
#> * Subject = 0 (population-level)
plot(pr)

When type = "re", the predicted values are still on the population-level. However, the random effect variances are taken into account, meaning that the intervals are actually prediction intervals and become larger. More technically speaking, type = "re" accounts for the uncertainty of the fixed effects conditional on the estimates of the random-effect variances and conditional modes (BLUPs).

The random-effect variance is the mean random-effect variance. Calculation is based on the proposal from Johnson et al. 2014, which is also implemented in functions like performance::r2() or insight::get_variance() to get r-squared values or random effect variances for mixed models with more complex random effects structures.

As can be seen, compared to the previous example with type = "fe", predicted values are identical (both on the population-level). However, standard errors, and thus the resulting confidence (or prediction) intervals are much larger .

pr <- ggpredict(m, "Days", type = "re")
pr
#> 
#> # Predicted values of Reaction
#> # x = Days
#> 
#>  x predicted std.error conf.low conf.high
#>  0   251.405    41.771  169.536   333.274
#>  1   261.872    41.765  180.015   343.730
#>  2   272.340    41.816  190.382   354.297
#>  3   282.807    41.924  200.638   364.976
#>  5   303.742    42.309  220.818   386.665
#>  6   314.209    42.584  230.745   397.673
#>  7   324.676    42.914  240.567   408.786
#>  9   345.611    43.730  259.902   431.319
#> 
#> Adjusted for:
#> * Subject = 0 (population-level)
plot(pr)

The reason why both type = "fe" and type = "re" return predictions at population-level is because ggpredict() returns predicted values of the response at specific levels of given model predictors, which are defined in the data frame that is passed to the newdata-argument (of predict()). The data frame requires data from all model terms, including random effect terms. This again requires to choose certain levels or values also for each random effect term, or to set those terms to zero or NA (for population-level). Since there is no general rule, which level(s) of random effect terms to choose in order to represent the random effects structure in the data, using the population-level seems the most clear and consistent approach.

To get predicted values for a specific level of the random effect term, simply define this level in the condition-argument.

ggpredict(m, "Days", type = "re", condition = c(Subject = 330))
#> 
#> # Predicted values of Reaction
#> # x = Days
#> 
#>  x predicted std.error conf.low conf.high
#>  0   275.094    41.771  193.225   356.963
#>  1   280.747    41.765  198.890   362.604
#>  2   286.400    41.816  204.443   368.358
#>  3   292.053    41.924  209.884   374.222
#>  5   303.360    42.309  220.436   386.284
#>  6   309.013    42.584  225.549   392.477
#>  7   314.666    42.914  230.556   398.776
#>  9   325.972    43.730  240.264   411.681

Finally, it is possible to obtain predicted values by simulating from the model, where predictions are based on simulate().

ggpredict(m, "Days", type = "sim")
#> 
#> # Predicted values of Reaction
#> # x = Days
#> 
#>  x predicted conf.low conf.high
#>  0   251.671  201.295   302.018
#>  1   261.947  212.139   311.333
#>  2   272.049  222.235   321.551
#>  3   282.964  232.901   333.247
#>  5   303.797  254.226   353.847
#>  6   314.486  264.274   364.580
#>  7   324.577  274.010   374.667
#>  9   345.632  295.128   395.344
#> 
#> Adjusted for:
#> * Subject = 0 (population-level)

Marginal effects for zero-inflated mixed models

For zero-inflated mixed effects models, typically fitted with the glmmTMB or GLMMadaptive packages, predicted values can be conditioned on

library(glmmTMB)
data(Salamanders)
m <- glmmTMB(
  count ~ spp + mined + (1 | site), 
  ziformula = ~ spp + mined, 
  family = truncated_poisson, 
  data = Salamanders
)

Similar to mixed models without zero-inflation component, type = "fe" and type = "re" for glmmTMB-models (with zero-inflation) both return predictions on the population-level, where the latter option accounts for the uncertainty of the random effects. In short, predict(..., type = "link") is called (however, predicted values are back-transformed to the response scale).

ggpredict(m, "spp")
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#>  x predicted std.error conf.low conf.high
#>  1     0.935     0.206    0.624     1.400
#>  2     0.555     0.308    0.304     1.015
#>  3     1.171     0.192    0.804     1.704
#>  4     0.769     0.241    0.480     1.233
#>  5     1.786     0.182    1.250     2.550
#>  6     1.713     0.182    1.200     2.445
#>  7     0.979     0.196    0.667     1.437
#> 
#> Adjusted for:
#> * mined = yes
#> *  site = NA (population-level)

ggpredict(m, "spp", type = "re")
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#>  x predicted std.error conf.low conf.high
#>  1     0.935     0.309    0.510     1.714
#>  2     0.555     0.384    0.261     1.180
#>  3     1.171     0.300    0.650     2.107
#>  4     0.769     0.333    0.400     1.478
#>  5     1.786     0.294    1.004     3.175
#>  6     1.713     0.294    0.964     3.045
#>  7     0.979     0.303    0.541     1.772
#> 
#> Adjusted for:
#> * mined = yes
#> *  site = NA (population-level)

For type = "fe.zi", the predicted response value is the expected value mu*(1-p) without conditioning on random effects. Since the zero inflation and the conditional model are working in “opposite directions”, a higher expected value for the zero inflation means a lower response, but a higher value for the conditional model means a higher response. While it is possible to calculate predicted values with predict(..., type = "response"), standard errors and confidence intervals can not be derived directly from the predict()-function. Thus, confidence intervals for type = "fe.zi" are based on quantiles of simulated draws from a multivariate normal distribution (see also Brooks et al. 2017, pp.391-392 for details).

ggpredict(m, "spp", type = "fe.zi")
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#>  x predicted std.error conf.low conf.high
#>  1     0.138     0.046    0.050     0.226
#>  2     0.017     0.009    0.000     0.035
#>  3     0.245     0.077    0.097     0.393
#>  4     0.042     0.018    0.006     0.078
#>  5     0.374     0.108    0.170     0.578
#>  6     0.433     0.121    0.199     0.667
#>  7     0.205     0.062    0.085     0.325
#> 
#> Adjusted for:
#> * mined = yes
#> *  site = NA (population-level)

For type = "re.zi", the predicted response value is the expected value mu*(1-p), accounting for the random-effect variances. Intervals are calculated in the same way as for type = "fe.zi", except that the mean random effect variance is considered and thus prediction intervals rather than confidence intervals are returned.

ggpredict(m, "spp", type = "re.zi")
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#>  x predicted std.error conf.low conf.high
#>  1     0.138     0.235    0.032     0.355
#>  2     0.017     0.231    0.000     0.054
#>  3     0.245     0.243    0.064     0.613
#>  4     0.042     0.231    0.005     0.118
#>  5     0.374     0.254    0.107     0.911
#>  6     0.433     0.259    0.134     1.031
#>  7     0.205     0.239    0.052     0.516
#> 
#> Adjusted for:
#> * mined = yes
#> *  site = NA (population-level)

Finally, it is possible to obtain predicted values by simulating from the model, where predictions are based on simulate() (see Brooks et al. 2017, pp.392-393 for details). To achieve this, use type = "sim".

ggpredict(m, "spp", type = "sim")
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#>  x predicted std.error conf.low conf.high
#>  1     1.094     1.281        0     4.120
#>  2     0.295     0.669        0     2.305
#>  3     1.525     1.554        0     5.316
#>  4     0.538     0.956        0     3.098
#>  5     2.225     2.138        0     7.197
#>  6     2.288     2.072        0     7.087
#>  7     1.326     1.373        0     4.666
#> 
#> Adjusted for:
#> * mined = yes
#> *  site = NA (population-level)

Marginal effects for each level of random effects

Marginal effects can also be calculated for each group level in mixed models. Simply add the name of the related random effects term to the terms-argument, and set type = "re".

In the following example, we fit a linear mixed model and first simply plot the marginal effetcs, not conditioned on random-effect variances.

library(sjlabelled)
data(efc)
efc$e15relat <- as_label(efc$e15relat)
m <- lmer(neg_c_7 ~ c12hour + c160age + c161sex + (1 | e15relat), data = efc)
me <- ggpredict(m, terms = "c12hour")
plot(me)

Changing the type to type = "re" still returns population-level predictions by default. Recall that the major difference between type = "fe" and type = "re" is the uncertainty in the variance parameters. This leads to larger confidence intervals (i.e. prediction intervals) for marginal effects with type = "re".

me <- ggpredict(m, terms = "c12hour", type = "re")
plot(me)

To compute marginal effects for each grouping level, add the related random term to the terms-argument. In this case, confidence intervals are not calculated, but marginal effects are conditioned on each group level of the random effects.

me <- ggpredict(m, terms = c("c12hour", "e15relat"), type = "re")
plot(me)

Marginal effects, conditioned on random effects, can also be calculated for specific levels only. Add the related values into brackets after the variable name in the terms-argument.

me <- ggpredict(m, terms = c("c12hour", "e15relat [child,sibling]"), type = "re")
plot(me)

The most complex plot in this scenario would be a term (c12hour) at certain values of two other terms (c161sex, c160age) for specific levels of random effects (e15relat), so we have four variables in the terms-argument.

me <- ggpredict(
  m, 
  terms = c("c12hour", "c161sex", "c160age", "e15relat [child,sibling]"), 
  type = "re"
)
plot(me)

If the group factor has too many levels, you can also take a random sample of all possible levels and plot the marginal effects for this subsample of group levels. To do this, use term = "<groupfactor> [sample=n]".

m <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
me <- ggpredict(m, terms = c("Days", "Subject [sample=7]"), type = "re")
plot(me)

You can also add the observed data points for each group using rawdata = TRUE.

plot(me, rawdata = TRUE)

References

Brooks ME, Kristensen K, Benthem KJ van, Magnusson A, Berg CW, Nielsen A, et al. glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal. 2017;9: 378–400.

Johnson PC, O’Hara RB. 2014. Extension of Nakagawa & Schielzeth’s R2GLMM to random slopes models. Methods Ecol Evol, 5: 944-946. (doi: 10.1111/2041-210X.12225)