################################################
# 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 = ")")