## Abstract

The coefficient of determination *R*^{2} quantifies the amount of variance explained by regression coefficients in a linear model. It can be seen as the fixed-effects complement to the repeatability *R* (intra-class correlation) for the variance explained by random effects and thus as a tool for variance decomposition. The *R*^{2} of a model can be further partitioned into the variance explained by a particular predictor or a combination of predictors using semi-partial (part) *R*^{2} and structure coefficients, but this is rarely done due to a lack of software implementing these statistics. Here, we introduce `partR2`, an R package that quantifies part *R*^{2} for fixed effect predictors based on (generalized) linear mixed-effect model fits. The package iteratively removes predictors of interest and monitors the change in *R*^{2} as a measure of the amount of variance explained uniquely by a particular predictor or a set of predictors. `partR2` also estimates structure coefficients as the correlation between a predictor and fitted values, which provide an estimate of the total contribution of a fixed effect to the overall prediction, independent of other predictors. Structure coefficients are converted to the total variance explained by a predictor, termed ‘inclusive’ *R*^{2}, as the square of the structure coefficients times total *R*^{2}. Furthermore, the package reports beta weights (standardized regression coefficients). Finally, `partR2` implements parametric bootstrapping to quantify confidence intervals for each estimate. We illustrate the use of `partR2` with real example datasets for Gaussian and binomials GLMMs and discuss interactions, which pose a specific challenge for partitioning the explained variance among predictors.

## Introduction

Coefficients of determination *R*^{2} are of interest in the study of ecology and evolution, because they quantify the amount of variation explained by a linear model (Edwards et al., 2008). By doing so, they go beyond significance testing in putting effects in perspective of the phenotypic variance. *R*^{2} is expressed as a proportion of the total variance in the response, which represents a biologically relevant quantity if the total variation is representative for the total population (de Villemereuil et al., 2018). The total coefficient of determination quantifies the variance explained by all fixed effects together (marginal *R*^{2} *sensu* Nakagawa & Schielzeth (2013), also known as the total correlation coefficient, Watanabe (1960)).

However, it is often of interest to attribute explained variation to individual predictors. Semi-partial coefficients of determination, also known as part *R*^{2}, decompose the variance of *R*^{2} into components uniquely explained by individual predictors (Jaeger et al., 2017) or sets of predictors (Figure 1). The set of all predictors in the model yields the total proportion of variance explained by the fixed part of the model (total *R*^{2}). With correlations among predictors, it often happens that predictors in univariate regressions explain a large share of the variance, but do not show large part *R*^{2} if other correlated predictors are included in the model.

Structure coefficients provide a valuable addition to part *R*^{2} in the decomposition of the phenotypic variance (Nimon et al., 2008). Structure coefficients quantify the correlation between individual predictors and the linear predictor (**η** = **Xβ**, where **X** is the design matrix for fixed effects and **β** is a vector of regression coefficients). Predictors that correlate well with a response, but are fitted with collinear predictors may show large structure coefficients as they are correlated to the predicted response, but low part *R*^{2} as other predictors explain part of the same variance. Structure coefficients range from −1 to 1 with their absolute value expressing the correlation relative to a perfect correlation if a single predictor explains as much as the total fixed part of the model.

Structure coefficients are correlations and since the square of a correlation yields the variance explained, we can use structure coefficients to estimate the total variance explained by a predictor (Nimon et al., 2008). We call this the inclusive *R*^{2} of a predictor and calculate it as the squared structure coefficient, i.e. its contribution to the linear predictor independent of other predictors (Nimon et al., 2008) times the proportion of variance explained by the linear predictor (which is the ‘total’ marginal *R*^{2} of the model). As far as we are aware, inclusive *R*^{2} has not been implemented before, but it provides valuable insights into the structure of the variance explained (Figure 1).

Here, we introduce `partR2`, a versatile package for estimating part *R*^{2}, inclusive *R*^{2}, structure coefficients and beta weights from mixed-effects models. The analysis is similar to, but not identical, to commonality analysis (Ray-Mukherjee et al., 2014; Seibold & McPhee, 1979; Zientek & Thompson, 2006, Figure 1). We illustrate how to use `partR2` with real example datasets for Gaussian and binomial GLMMs, discuss how to estimate part R^{2} in the presence of interactions and discuss some challenges and limitations.

### Other implementations in R packages

There are a few R packages that calculate part *R*^{2} for linear models (lm), for example `rockchalk∷getDeltaRsquare` (Johnson & Grothendieck, 2019). Other packages calculate partial *R*^{2} (not part *R*^{2}) such as `asbio∷partial.R2` (Aho, 2020) and `rr2∷R2` (Ives & Li, 2018) for lms and `rsq∷rsq.partial` (Zhang, 2020) for linear models and generalized linear models (glm). Note that partial *R*^{2} is different from part (semi-partial) *R*^{2} (partial *R*^{2} > part *R*^{2}), since it represents the unique variance explained by a particular predictor but after removing (‘partialling out’) the variance explained by the other predictors (Yeatts et al., 2017, Figure 1). The package yhat features a suite of functions for full communality analyses in glms (Nimon, Oswald & Roberts, 2020). None of these packages estimates part *R*^{2} for mixed-effects models that we focus on here.

Several packages estimate (marginal) *R*^{2} as the variance explained by all fixed effects in linear mixed-effects models. This includes `performance∷r2_nakagawa` (Lüdecke et al., 2020), `MuMIn∷r.squaredGLMM` (Bartoń, 2019), and `rptR∷rpt` (Stoffel, Nakagawa & Schielzeth, 2017). These packages do not allow to estimate part *R*^{2}. The only versatile package to estimate part *R*^{2} from linear mixed-models is `r2glmm` (Jaeger, 2017). The function `r2glmm∷r2beta` computes part *R*^{2} from lmer, lme and glmmPQL model fits (also for linear models lm and glm) based on Wald statistics. However, it does neither support `lme4∷glmer` for generalized linear model fits nor does it allow to estimate *R*^{2} for combinations of predictors. Furthermore, it does not estimate structure coefficients, inclusive *R*^{2} or part *R*^{2} for multilevel factors as a unit.

### 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*^{2} by iterative removal of 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*^{2} 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*^{2} *(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*^{2}, inclusive *R*^{2}, structure coefficients, beta weights (standardized regression slopes), bootstrapping iterations and some other information. An extended summary, that includes inclusive *R*^{2}, 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., in prep.).

*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*^{2} 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, 2019). 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”),
data=GuineaPigs, 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*^{2}, structure coefficients and beta weights (standardized slopes) (Figure 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*^{2} (with point estimates and confidence intervals), `res$SC` for structure coefficients, `res$IR2` for inclusive *R*^{2} and `res$BW` for beta weights.

### Dealing with interactions

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

The model output above shows the number of parameters fitted in each model (Figure 2, 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*^{2} 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, while 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*^{2} 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”)), data=GuineaPigs,
nboot=100)

This, however, is only one way of dealing with interactions (Option A in Figure 3). It represents the variance explained uniquely 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*^{2} 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 Figure 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”)), data=GuineaPigs, 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 Figure 3). This implies that part *R*^{2} 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”), data=GuineaPigs,
nboot=100)
mod2 <- lmer(Testo ~ Time + Rank + (1|MaleID), data=GuineaPigs)
part2 <- partR2(mod2, partvars = c(“Time”, “Rank”), data=GuineaPigs,
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 Figure 3.

### An example with proportion data

As an example for proportion data, we analyze a dataset on spatial variation in color morph ratios in a color-polymorphic species of grasshopper. Individuals of this species occur either in a green or a brown color variant and the dataset contains counts of brown and green individuals (separated for females and males) from 42 sites sampled in the field (Dieker et al., 2018). Site identity will be fitted as a random effect. As covariates the dataset contains a range of Bioclim variable that describe various aspects of ecologically relevant climatic conditions (see Karger *et al.* 2017). The aim is to identify the climatic conditions that favour one or the other colour variant.

We first load the grasshopper dataset. We standardise all Bioclim variables using the `scale` function and add an observation-level counter that will be used as an observation-level random effect (OLRE) to account for overdispersion (Harrison, 2014).
data(Grasshoppers)
for (i in which(substr(colnames(Grasshoppers),1,3)==“Bio”)){
Grasshoppers[,i] <- scale(Grasshoppers[,i])
}
Grasshoppers$OLRE <- 1:nrow(Grasshoppers)

We first fit a GLMM with binomial error structure and logit link using the `glmer` function from the `lme4` package (Bates et al., 2015). A previous analysis has shown that the first principle component of the Bioclim data explains a small, but significant part of variation in morph ratios (Dieker et al., 2018). For illustration, we use the four Bioclim variables that show a loading of more than 0.30 on the first principle component.
mod <- glmer(cbind(nGreen, nBrown) ~ Bio7 + Bio14 + Bio17 + Bio19 +
(1|SiteID) + (1|OLRE), data=Grasshoppers, family=“binomial”)
res <- partR2(mod, partvars=c(“Bio7”, “Bio14”, “Bio17”, “Bio19”),
data=Grasshoppers, max_level = 1, nboot=100)

The `summary` output informs us (at the bottom) that there have been warnings in the bootstrapping processes. This is not unusual since bootstrapping frequently generates data, for which one of the parameters is estimated at the boundary (in particular if one of the variance components is very small). The results can be visualised using the `forestplot` function (Figure 4). Plotting is based on `ggplot2` (Wickham, 2016), and multiple forest plots can easily be assembled using the `patchwork` package (Pedersen, 2020). Forest plots show the effect sizes graphically and can be set to either show part *R*^{2} when `type = “R2”` (the default), inclusive *R*^{2} when `type = “IR2”`, structure coefficients when `type = “SC”`, and beta weights (standardized model estimates) with `type = “BW”`.
p1 <- forestplot(res, type = “R2”)
p2 <- forestplot(res, type = “IR2”)
p3 <- forestplot(res, type = “SC”)
p4 <- forestplot(res, type = “BW”)
library(patchwork)
(p1 + p2) / (p3 + p4) + plot_annotation(tag_levels = “A”, tag_prefix =
“(“, tag_suffix = ”)”)

A comparison of part *R*^{2}, inclusive *R*^{2}, structure coefficients beta weights shows the different insights that can be gained from these different summaries of the model fit (Figure 3). In this case, three of the Bioclim variables (*Bio14*, *Bio17*, *Bio19*) are highly positively correlated (r ≥ 0.93), while a fourth one (*Bio7*) is moderately negatively correlated to all three of them (r ≤ −0.63). Part *R*^{2} are thus low, because none of the parameters uniquely explains a large share of the variance. *Bio17* seems to be the best predictor of morph ratios, with the largest (negative) beta weight, largest part *R*^{2}, largest structure coefficients and largest inclusive *R*^{2}. Beta weights for the two positively correlated (but slightly weaker) predictors, *Bio14* and *Bio19*, switch sign as is not unusual for collinear predictors.

This means that after accounting for the effect of *Bio17*, they contribute positively to prediction. However, structure coefficients show that both variables load negatively on the linear predictor, as does *Bio17*.

### 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, 2007). 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*^{2}, 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*^{2} (see also Rights & Sterba, 2019). By default, partR2 sets negative *R*^{2} values to 0, but this can be changed by setting `allow_neg_r2` to `TRUE`. It also happens that inclusive *R*^{2} is estimated slightly lower than part *R*^{2} 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*^{2} (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 development version of `partR2` can be obtained from GitHub (https://github.com/mastoffel/partR2) and includes the data used in the examples. The package will be uploaded to CRAN after peer-review.

## Author contributions

SN and HS conceived the idea. HS and MAS wrote the manuscript with contributions from SN. MAS wrote the package and the vignette, supervised by HS. All authors contributed critically to the manuscript, package and vignette and gave final approval for publication.

## Acknowledgements

MAS and HS were supported by the German Research Foundation (DFG) as part of the SFB TRR 212 (NC³) (funding INST 215/543-1, 396782608). SN was supported by the ARC Discovery Project grant (DP180100818). We thank Alexandra Mutwill for providing the guinea pig data.