partR2: partitioning R² in generalized linear mixed models

Part R²

A Gaussian mixed-effects model can be written as: (1)${\displaystyle \mathbf{y}=\mathbf{X}\beta +\sum {\alpha }_{\mathrm{k}}+\epsilon }$ ${\displaystyle {\alpha }_k\sim N\left(0,{\sigma }_{{\alpha }_k}^2\right)}$ ${\displaystyle \epsilon \sim N\left(0,{\sigma }_{\epsilon }^2\right)}$

Where y is a vector of response values (outcomes), X is the design matrix of fixed effects, β is a vector of regression coefficients, ∑α_k is the random part of the model that might contain multiple random effects and ε is a vector of residual deviations. The linear predictor η represents the vector of predicted values from the fixed part of the model as η = Xβ. Note that we dealing with estimates of regression coefficients and variance components throughout (hence all β should be read as $\widehat{\beta }$).

Since we are interested in the proportion of the phenotypic variance explained, we symbolize variance components by upper case Y and index by the source of variance (as in Fig. 1). While variances are frequently represented as V with the source of variance as an index, this leads to ambiguity for V_X which might represent variance in y explained by x or variance in x itself, which is why we use this alternative notation. The total variance in the response is Y_Total = var(y) and is estimated from the raw data or from the model (see below). The variance of the residuals is estimated by the model as Y_R = var(ε). The variance of the (sum of) random effects is estimated by the model as Y_RE = var(∑α_k) and the variance explained by fixed effects can be estimated as the variance in the linear predictor Y_X = var(Xβ).

The coefficient of determination R² estimates the proportion of variance in the response that is explained by fixed effects. The coefficient of determination R² for the total fixed part of the model is thus: (2)${\displaystyle R_{{}^X}^2=\frac{Y_X}{Y_X+Y_{RE}+Y_R}=\frac{Y_X}{Y_{Total}}}$

Note that the sum of the components in the denominator might deviate numerically from the total outcome variance in the raw data. However, conceptually they are the same in that they represent the population-level outcome variance. The variance in the outcome is an estimate from the specific sample, while the sum of components of the mixed model represents a population-level estimate given the data and the model.

A reduced model with a (set of) fixed effect predictors X^∗ removed but the same random effect structure can be fitted as (now using the tilde to highlight the differences from Eq. (1)): (3)${\displaystyle \mathbf{y}=\tilde{\mathbf{X}}\tilde{\beta }+\sum {\tilde{\alpha }}_k+\tilde{\epsilon }}$ ${\displaystyle {\tilde{\alpha }}_k\sim N\left(0,{\sigma }_{{\tilde{\alpha }}_k}^2\right)}$ ${\displaystyle \tilde{\epsilon }\sim N\left(0,{\sigma }_{\tilde{\epsilon }}^2\right)}$ with the variance in the linear predictor of the reduced model being $Y_{\tilde{X}}=var\left(X\beta \right)$.

The variance uniquely explained by X^∗ is then the difference between the variance explained by fixed effects in the full and the reduced model $Y_{X^{\ast }}=Y_X-Y_{\tilde{X}}$. Part R² sets this variance in proportion to the total outcome variance: (4)${\displaystyle R_{X^{\ast }}^2=\frac{Y_X-Y_{\tilde{X}}}{Y_X+Y_{RE}+Y_R}=\frac{Y_X-Y_{\tilde{X}}}{Y_{Total}}}$

The process of fitting a reduced model, estimation of Y_X^∗ and estimation of $R_{{}^{X^{\ast }}}^2$ can be repeated for all predictors and combinations of predictors. At the limit for a model with all fixed effects removed, $R_{X\ast }^2=R_X^2$.

Inclusive R²

Structure coefficients are the Pearson correlations between a particular predictor of interest x^∗ and the linear predictor η. Note that we now use a lower case x^∗ to indicate that we are dealing with a single predictor. Structure coefficients are quantified from the full model as: (6)${\displaystyle S{C}_{x^{\ast }}=\mathrm{cor}\left(\eta ,x^{\ast }\right)}$

The squared correlation between two variables a and b gives the variance explained for these variables $\mathrm{cor}{\left(a,b\right)}^2=R_a^2$. The squared structure correlations thus quantify the proportion of variance in the linear predictor Y_X that is explained by a the predictor of interest x^∗. Since the proportion of outcome variance explained by the linear predictor in the full model is $R_{{}^X}^2$, the inclusive variance explained by predictor x^∗ is: (7)${\displaystyle I{R}_{x^{\ast }}^2=S{C}^2\cdot R_{{}^{X^{\ast }}}^2}$

Inclusive R² as we define it here, complements part R² by giving additional insights. While part R² quantifies the variance uniquely explained by a predictor (or set of predictors), inclusive R² quantifies the total proportion of variance explained in the model, both uniquely and jointly with other predictors. In the special case of a single predictor in a model SC_x^∗ = cor(η, x^∗) = 1, such that $I{R}_{x^{\ast }}^2=R_X^2$.

Part R² in non-Gaussian models

For Gaussian models there is a single residual error term ε with variance Y_R = var(ε). For non-Gaussian models, however, there is additional error that arises from the link function that translates latent-level predictions to observed outcomes. This variance can be approximated for a variety of link functions and error distributions (Nakagawa & Schielzeth, 2010; Nakagawa, Johnson & Schielzeth, 2017). Our R package currently implements distribution-specific variances for Poisson models with log and square root link functions and binomial models with logit and probit link functions. For Poisson models and non-binary binomial models (proportion models), partR2 also fits an observational level random effect (if none is fitted already) to estimate variance due to overdispersion (Harrison, 2014). Both the overdispersion variance, now denoted Y_R and the distribution-specific variance Y_D are included in the denominator of the part R2 calculation: (8)${\displaystyle R_{X\ast }^2=\frac{Y_X-Y_{\tilde{X}}}{Y_X+Y_{RE}+Y_D+Y_R}}$

Notably, there are other estimation methods for R² for non-Gaussian models or GLMM (Jaeger et al., 2017; Piepho, 2019). Currently, partR2 only implements the method based on Nakagawa & Schielzeth (2013) and Nakagawa, Johnson & Schielzeth (2017).

Features of partR2

partR2 takes a fitted (generalized) linear mixed-model (GLMM), from the popular mixed model package lme4 (Bates et al., 2015) and estimates part R² by iteratively removing fixed effects (Nimon et al., 2008). The specific fixed effects of interest are specified by the partvars and/or by the partbatch argument. The package estimates part R² for all predictors specified in partvars individually and in all possible combinations (the maximum level of combinations can be set by the max_level argument). A custom specification of fixed effects of interest saves computation time as compared to an all-subset specification and is therefore required in partR2.

The central function partR2 will work for Gaussian, Poisson and binomial GLMMs. Since the model fit is done externally, there is no need to supply a family argument. For non-Gaussian GLMMs, the package estimates link-scale R² (sensu Nakagawa & Schielzeth, 2013). We implement parametric bootstrapping to quantify sampling variance and thus uncertainty in the estimates. Parametric bootstrapping works through repeated model fitting on simulated data based on fitted values (Faraway, 2015). The number of bootstrap iterations is controlled by the nboot argument. We recommend a low number of nboot for testing purposes and a large number (e.g., nboot = 1000) for the final analysis.

The package returns an object of class partR2 that contains elements for part R², inclusive R², structure coefficients, beta weights (standardized regression slopes), bootstrapping iterations and some other information. An extended summary, that includes inclusive R², structure coefficients and beta weights can be viewed using the summary function. The forestplot function shows a graphical representation of the variance explained by individual predictors and sets of predictors along with their bootstrapping uncertainties. All computations can be parallelized across many cores based on the future and furrr packages (Vaughan & Dancho, 2018; Bengtsson, 2020). An extended vignette with details on the complete functionality accompanies the package.

Example with Gaussian data

We use an example dataset with hormone data collected from a population of captive guinea pigs to illustrate the features of partR2. The dataset contains testosterone measurements of 31 male guinea pigs, each measured at 5 time points (age between 120 and 240 days at 30-day intervals). We analyze log-transformed testosterone titers and fit male identity as a random effect. As covariates the dataset contains the time point of measurement and a rank index derived from behavioral observations around the time of measurement (Mutwill et al., 2021).

Rank and Time are correlated in the dataset (r = 0.40), since young individuals are typically low rank, while older individuals tend to hold a high rank. Time might be fitted as a continuous predictor or as a factor with five levels. Here we present the version of a factorial predictor to illustrate the estimation of part R² for interactions terms. Hence, an interaction between Time and Rank will also be fitted.

First, the package needs to be loaded (after successful installation) in an R session (R Core Team, 2021). The package comes with the guinea pig dataset that also needs to be loaded using the data function.

library(partR2)

data(GuineaPigs)

A single record contains missing values for testosterone measurements. Missing records can be problematic to handle in partR2 and are better removed prior to the analysis. We also log-transform the response and convert Time to a factor and filter for the first three time points to simplify the output.

GuineaPigs <- subset(GuineaPigs, !is.na(Testo) & !is.na(Rank) & (Time %in% c(1,3,5)))

GuineaPigs$TestoTrans <- log(GuineaPigs$Testo)

GuineaPigs$Time <- factor(GuineaPigs$Time)

We then fit a linear mixed effects model using lmer from the lme4 package (Bates et al., 2015). Further exploration of the data and model checks are omitted here for simplicity, but are advisable in real data analysis.

library(lme4)

mod <- lmer(TestoTrans ∼Rank * Time + (1—MaleID), data=GuineaPigs)

The partR2 analysis takes the lmer model fit (an merMod object) and a character vector partvars indicating the fixed effects to be evaluated. Interactions are specified with the colon syntax (see the package’s vignette for further details).

res <- partR2(mod, partvars = c(”Rank”, ”Time”, ”Rank:Time”), nboot=100)

The function returns a partR2 object. The print function reports the part coefficients of determination and a more extensive summary can be viewed with the summary function which also shows inclusive R², structure coefficients and beta weights (standardized slopes) (Fig. 2).

print(res)

summary(res, round_to = 2)

The variances appear largely additive, since combinations of predictors explain about the sum of the variance explained by individual predictors. The main components of the partR2 object can be accessed for further processing as res$R2 for part R² (with point estimates and confidence intervals), res$SC for structure coefficients, res$IR2 for inclusive R² and res$BW for beta weights.

Dealing with interactions

Models with interactions are problematic, because the variance explained by a main factor can be estimated in multiple ways (Fig. 3) and because of the internal parametrization of the model matrix.

The model output in Fig. 2 shows the number of parameters fitted in each model (each row in the R2 part refers to a reduced model). In the print and summary output this is visible as a column labelled ‘ndf’. A close inspection shows that the removal of rank did not change the number of parameters (6 for the full model, 6 for the model excluding rank). This is because the model matrix is reparametrized in the reduced model and lmer will fit three terms for the interaction (here Time1:Rank, Time3:Rank, Time5:Rank) rather than just two for the interaction in the full model. Dummy coding of the factor can be usefully combined with centering of dummy coded variables (Schielzeth, 2010) and gives more control over this re-parametrisation. It allows for example to estimate the part R² for the average effect of Rank by constraining the average Rank effect to zero, so that only the two contrasts are fitted (here Time3:Rank, Time5:Rank):

GuineaPigs <- cbind(GuineaPigs, model.matrix(∼0 + Time, data=GuineaPigs))

GuineaPigs$Time3 <- GuineaPigs$Time3 - mean(GuineaPigs$Time3)

GuineaPigs$Time5 <- GuineaPigs$Time5 - mean(GuineaPigs$Time5)

The model can then be fitted with dummy predictors. Since the usual specification in partR2 via partvars would fit all possible combinations, including combinations of the different Time terms, such a run can take a long time. However we are mostly interested in fitting and removing all dummy predictors at a time. The package therefore features an additional argument partbatch to specify a list of character vectors containing the sets of predictors that should always be kept together. In the example, the list has two elements, a character vector for the dummy-coded main effects and a character vector for the interaction terms. The analysis yields part R² for two batches of predictors as well as Rank and their combinations.

mod <- lmer(TestoTrans ∼(Time3 + Time5) * Rank + (1—MaleID),

data=GuineaPigs)

batch <- c(”Time3”, ”Time5”)

partR2(mod, partvars=c(”Rank”), partbatch=list(Time=batch, “Time:Rank”= paste0(batch, ”:Rank”)), nboot=100)

This, however, is only one way of dealing with interactions (Option A in Fig. 3). It represents the variance uniquely explained by main effects even in the presence of an interaction. Since interactions are the products of main effects, interaction terms are typically correlated with main effects and the part R² calculated above might not represent a biologically relevant quantity. There are two alternative ways of how to deal with interactions. Both are possible in partR2, but since requirements differ between applications, we do not implement one with priority.

One way to think about variance explained by main effects and their interactions is to pool the variance explained by a main effect with the variance explained by interactions that the term is involved in (Option B in Fig. 3). In the guinea pig example, for instance, Rank might be considered important either as a main effect or in interaction with time and we might want to estimate the total effect of Rank. This can be done for the guinea pig dataset by using partbatch:

mod <- lmer(Testo ∼Time * Rank + (1—MaleID), data=GuineaPigs)

partR2(mod, partbatch = list(Time=c(”Time”, ”Time:Rank”), Rank=c(”Rank”, ”Time:Rank”)), nboot=100)

A third, which we think usually preferable option is to prioritize main effects by assigning the proportion of variance that is explained by a main effect together with the variance jointly explained with its interaction to the main effect (Option C in Fig. 3). This implies that part R² for a main effect is estimated when its own interaction is excluded from the model (mod1 and part1 below). The variance explained by the interaction is then estimated in a separate model (mod2 and part2 below). We have implemented a helper function mergeR2 that allows to merge two partR2 runs.

mod1 <- lmer(Testo ∼Time * Rank + (1—MaleID), data=GuineaPigs)

part1 <- partR2(mod1, partvars = c(”Time:Rank”), nboot=100)

mod2 <- lmer(Testo ∼Time + Rank + (1—MaleID), data=GuineaPigs)

part2 <- partR2(mod2, partvars = c(”Time”, ”Rank”), nboot=100)

mergeR2(part1, part2)

All these results can be viewed by print, summary and plotted by forestplot. It is important to bear in mind the differences in the interpretation as illustrated in Fig. 3.

Challenges

Using transformation or functions in the formula argument can lead to issues with matching the terms of the model with the partvars argument of partR2. It is therefore important that the names in partvars match exactly the terms in the merMod object. However, any complications are easily circumvented by implementing the transformations before fitting the model and storing them in the data frame used in the analysis. It is also worth to be aware that unusual names may cause complications and renaming can offer an easy solution.

We have repeatedly seen model outputs where the point estimate does not fall within the confidence interval. This might seem like in the bug in the package, but in our experience usually indicates issues with the data and/or the model. In fact, parametric bootstrapping can be seen as a limited form of posterior predictive model checks (Gelman & Hill, 2006). If generating new data from the fitted model (as done with parametric bootstrapping) results in data that are dissimilar to the original data, then the model is probably not a good fit to the data.

Bootstrap iterations can sometimes yield slightly negative estimates of part R², in particular if the variance explained by a predictor is low. These negative estimates happen in mixed-effects models, because estimates of random-effect variance might change when a predictor is removed and this can lead to a slight decrease in the residual variance, and hence a proportional increase in R² (Rights & Sterba, 2019). By default, partR2 sets negative R² values to 0, but this can be changed by setting allow_neg_r2 to TRUE. It also happens that inclusive R² is estimated slightly lower than part R² when the contribution of a particular predictor is very large. We consider both cases as sampling error that should serve as a reminder that variance components are estimated with relatively large uncertainly and minor differences should not be over-interpreted.

A warning needs to be added for the estimation of R² (and, in fact, also repeatability R) from small datasets. In particular if the number of levels of random effect is low, variance components might be slightly overestimated (Xu, 2003). This issue applies similarly to the variance explained by fixed effects, in particular if the number of predictors is large relative to the number of data points.

Code and data availability

The current stable version of partR2 can be downloaded from CRAN (https://cran.r-project.org/web/packages/partR2/index.html) and the development version can be obtained from GitHub (https://github.com/mastoffel/partR2). The data used in the examples is part of the package.