Example Code for Calculating R2LMM(M) and R2LMM(C)

Supplementary Material

Example Code for Calculating R2LMM(M) and R2LMM(C)

setwd("Your working directory path") # This needs to be updated to whatever directory the data is located on your computer

#This code includes how to calculate marginal and conditional R2 values, as well as p-values for a linear mixed model.

NH4.DATA<-read.csv(file="NH4.DATA.csv", header=T) #This .csv files being read in needs to have Vf, Uvol, and all potential predictors for each chamber

#Note - before proceeding with the code you should check assumptions of different tests (code not included here)

#Load in packages necessary for the code:

library(arm)

library(lme4)

#First we will calculate the marginal and conditional Rsquared:

#We first write a function to calculate R2LMM

# RsquareGLME code from jonlefcheck.net blog post:

#

#last accessed in Oct 2014. The code was modified from the original code included

#in Nakagawa and Schielzeth 2013.

rsquared.glme=function(modlist) {

# Iterate over each model in the list

do.call(rbind,lapply(modlist,function(i) {

# For models fit using lm

if(class(i)=="lm") {

Rsquared.mat=data.frame(Class=class(i),Family="Gaussian",

Marginal=summary(i)$r.squared,Conditional=NA,AIC=AIC(i)) }

# For general linear models fit using lme4

else if(class(i)=="lmerMod" | class(i)=="merLmerTest") {

# Get variance of fixed effects by multiplying coefficients by design matrix

VarF=var(as.vector(fixef(i) %*% t(i@pp$X)))

# Get variance of random effects by extracting variance components

VarRand=colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))

# Get residual variance

VarResid=attr(VarCorr(i),"sc")^2

# Calculate marginal R-squared (fixed effects/total variance)

Rm=VarF/(VarF+VarRand+VarResid)

# Calculate conditional R-squared (fixed effects+random effects/total variance)

Rc=(VarF+VarRand)/(VarF+VarRand+VarResid)

# Bind R^2s into a matrix and return with AIC values

Rsquared.mat=data.frame(Class=class(i),Family="Gaussian",Marginal=Rm,Conditional=Rc,

AIC=AIC(update(i,REML=F))) }

#For generalized linear models (family=="binomial") fit using lme4

else if(class(i)=="glmerMod" & summary(i)$family=="binomial") {

# Get variance of fixed effects by multiplying coefficients by design matrix

VarF=var(as.vector(fixef(i) %*% t(i@pp$X)))

# Get variance of random effects by extracting variance components

VarRand=colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))

# Get residual variance

VarResid=attr(VarCorr(i),"sc")^2

# Calculate marginal R-squared

Rm=VarF/(VarF+VarRand+pi^2/3)

# Calculate conditional R-squared (fixed effects+random effects/total variance)

Rc=(VarF+VarRand)/(VarF+VarRand+pi^2/3)

# Bind R^2s into a matrix and return with AIC values

Rsquared.mat=data.frame(Class=class(i),Family=summary(i)$family,Marginal=Rm,Conditional=Rc,AIC=AIC(i)) }

#For generalized linear models (family=="poisson") fit using lme4

else if(class(i)=="glmerMod" & summary(i)$family=="poisson") { print("GLMM fit to Poisson dist not yet supported") }

# For model fit using nlme

else if(class(i)=="lme") {

# Get design matrix of fixed effects from model

Fmat=model.matrix(eval(i$call$fixed)[-2],i$data)

# Get variance of fixed effects by multiplying coefficients by design matrix

VarF=var(as.vector(fixef(i) %*% t(Fmat)))

# Get variance of random effects by extracting variance components

VarRand=sum(suppressWarnings(as.numeric(VarCorr(i)[rownames(VarCorr(i))!=

"Residual",1])),na.rm=T)

# Get residual variance

VarResid=as.numeric(VarCorr(i)[rownames(VarCorr(i))=="Residual",1])

# Calculate marginal R-squared (fixed effects/total variance)

Rm=VarF/(VarF+VarRand+VarResid)

# Calculate conditional R-squared (fixed effects+random effects/total variance)

Rc=(VarF+VarRand)/(VarF+VarRand+VarResid)

# Bind R^2s into a matrix and return with AIC values

Rsquared.mat=data.frame(Class=class(i),Marginal=Rm,Conditional=Rc,

AIC=AIC(update(i,method="ML")))

} else { print("Function requires models of class lm, lme, mer, or merMod")

} } ) ) }

#Next we create a null and alternative model

nh4.null.model<-lmer(LOG.NH4.VF~1+(1|SITE), data=NH4.DATA)

nh4.lmm<-lmer(LOG.NH4.VF~LOG.WS.AREA*WS+(1|SITE), data=NH4.DATA)

#And we find the marginal and conditional R2 values for each model:

rsquared.glme(list(nh4.null.model, nh4.lmm))

#Next we use lmerTest to get P-values

#Note that the lmerTest package makes calculating LMMR2s impossible, therefore

#p-values must be calculated after R2s.

library(lmerTest)

#Re-run your LMM

nh4.lmm<-lmer(LOG.NH4.VF~(LOG.WS.AREA+WS)^2+(1|SITE), data=NH4.DATA)

#Get the p-values for the model

anova(nh4.lmm)