################################################
# Installation
################################################
# the development version of partR2 can be downloaded from GitHub
# using the "remotes" package. After peer review, the package will be
# uploaded to CRAN
# install.packages("remotes")
# remotes::install_github("mastoffel/partR2", build_vignettes = TRUE)
library(partR2)
# vignette
browseVignettes(package = "partR2")
################################################
# Example with Gaussian data: Guinea Pig dataset
################################################
# load dataset
data(GuineaPigs)
# subset non-missing values for testosterone and rank and subset
# only three time points for simpler outputs
GuineaPigs <- subset(GuineaPigs, !is.na(Testo) & !is.na(Rank) & (Time %in% c(1,3,5)))
# log transform testosterone measurements
GuineaPigs$TestoTrans <- log(GuineaPigs$Testo)
# treat time as factor
GuineaPigs$Time <- factor(GuineaPigs$Time)
# run linear mixed model with Rank and Time main effects and interaction
# using *
library(lme4)
mod <- lmer(TestoTrans ~ Rank * Time + (1|MaleID), data=GuineaPigs)
# set a seed to get the same output as in the paper
set.seed(123)
# partR2 analysis. We partition the R2 for Rank and Time main effects and
# their interaction and use 100 parametric bootstraps to estimate confidence
# intervals
res <- partR2(mod, partvars = c("Rank", "Time", "Rank:Time"), nboot=100)
# print and summary give short and more extensive results
print(res)
summary(res, round_to = 2)
################################################
# Dealing with interactions
################################################
# Transforming the factor Time into a matrix of dummy variables (0/1) for each
# time point using model.matrix() and attaching these dummy variables
# to the GuineaPigs data.frame
GuineaPigs <- cbind(GuineaPigs, model.matrix(~ 0 + Time, data=GuineaPigs))
# centering the time dummy variables which we will use for modeling
GuineaPigs$Time3 <- GuineaPigs$Time3 - mean(GuineaPigs$Time3)
GuineaPigs$Time5 <- GuineaPigs$Time5 - mean(GuineaPigs$Time5)
# (1) first way to estimate part R2 for interactions (Figure 3A in manuscript)
# note that in contrast to the generic syntax in Figure 3A, this uses
# the partbatch argument to remove all "Time" terms at once, which can be
# useful for models with many variables and models with dummy variables
# to prevent fitting to many models which would take a long time.
# this gets the variance uniquely explained by the main effects and by
# the interaction without their 'shared' explained variance.
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)
# (2) second way to estimate part R2 for interactions (Figure 3B in manuscript)
# We create a list for partbatch. The first element, called Time will be used
# to estimate the combined part R2 of the Time main effect and the Time:Rank
# interaction. Similarly, the second element, called Rank will estimate the
# combined part R2 of the Rank main effect the Time:Rank interaction.
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)
# (3) third way to estimate part R2 for interactions (Figure 3C in manuscript)
# usually the preferred way
# This assigns the jointly explained variance from a main effect and its
# interaction to the main effect part R2. Here, we need to fit two models
# and combine them:
# model one contains main effects and their interaction
mod1 <- lmer(Testo ~ Time * Rank + (1|MaleID), data=GuineaPigs)
# Create partR2 object by estimating partR2 for the Time:Rank interaction
part1 <- partR2(mod1, partvars = c("Time:Rank"), nboot=100)
# model two contains only the main effects
mod2 <- lmer(Testo ~ Time + Rank + (1|MaleID), data=GuineaPigs)
# Create partR2 object by estimating partR2 for only the main effects
part2 <- partR2(mod2, partvars = c("Time", "Rank"), nboot=100)
# now use mergeR2 to combine both into a new partR2 object
mergeR2(part1, part2)
################################################
# Example with Proportion data: Grasshopper dataset
################################################
data(Grasshoppers)
# scale() function standardises all Bioclim variables
# Bioclim variables describe ecologically relevant climatic conditions
# which might impact color morph ratios
for (i in which(substr(colnames(Grasshoppers),1,3)=="Bio")){
Grasshoppers[,i] <- scale(Grasshoppers[,i])
}
# Adding an observational level random effect OLRE to model overdispersion
Grasshoppers$OLRE <- 1:nrow(Grasshoppers)
# model proportion data using binomial GLMM
# response is the color morph ratio and fixed effects are
# the Bioclim variables. Random effects are the sampling site ID and
# the OLRE
mod <- glmer(cbind(nGreen, nBrown) ~ Bio7 + Bio14 + Bio17 + Bio19 + (1|SiteID) + (1|OLRE),
data=Grasshoppers, family="binomial")
# partR2 analysis. Here we use max_level = 1, which gives us partR2
# only for each of the predictors but not for their combinations.
res <- partR2(mod, partvars=c("Bio7", "Bio14", "Bio17", "Bio19"),
nboot=100, max_level = 1)
# plotting all partR2 outputs as forestplot
p1 <- forestplot(res, type = "R2")
p2 <- forestplot(res, type = "IR2")
p3 <- forestplot(res, type = "SC")
p4 <- forestplot(res, type = "BW")
# Each forestplot is a ggplot. The patchwork package can be used to easily
# combine ggplots and annotate them.
library(patchwork)
(p1 + p2) / (p3 + p4) + plot_annotation(tag_levels = "A", tag_prefix = "(", tag_suffix = ")")