########################################## ############ step 1. get set up to work ######## ########################################## library(lme4) library(afex) library(ENmisc) library(arm) ########################################## ############ step 2. import datafile ######## ########################################## sis <-read.csv ("sister.summary.table.csv") head(sis) names(sis) ########################################## ############ step 3. drop unwanted levels ######## ########################################## sis.sub=subset(sis,mate.diff=="s_o") levels(sis.sub$mate.diff) sis.sub$mate.diff <- factor(sis.sub$mate.diff) levels(sis.sub$mate.diff) nrow(sis.sub) ########################################## ############ step 4. restructure data ######## ########################################## test=subset(sis.sub, mate.diff=="s_o") test$mate.diff <- factor(test$mate.diff) levels(test$mate.diff) test$range.diff=test$A.range.res1 - test$B.range.res1 self=subset(test, select=c(A.range.res1,A.range.res0.5,A.range.res0.1,A.range.res0.05,pp,fig.weight,genus,dist,ba,lat.med.spA,lat.min.spA,lat.max.spA,lat.mean.spA)) out=subset(test, select=c(B.range.res1,B.range.res0.5,B.range.res0.1,B.range.res0.05,pp,fig.weight,genus,dist,ba,lat.med.spB,lat.min.spB,lat.max.spB,lat.mean.spB)) self$mate.sys="selfer" out$mate.sys="outcrosser" colnames(self)=c("range.size.res1","range.size.res0.5","range.size.res0.1","range.size.res0.05","pp","fig.weight","genus","dist","ba","lat.med","lat.min","lat.max","lat.mean","mate.sys") colnames(out)=c("range.size.res1","range.size.res0.5","range.size.res0.1","range.size.res0.05","pp","fig.weight","genus","dist","ba","lat.med","lat.min","lat.max","lat.mean","mate.sys") tmp3=rbind(self,out) tmp3$mate.sys = factor(tmp3$mate.sys,c("selfer","outcrosser")) ########################################## ############ step 5. get set up for making figures ######## ########################################## library(RColorBrewer) brewer.pal(n=20, name="Set1") brewer.pal(n=20, name="Set2") brewer.pal(n=20, name="Set3") par(mfrow=c(1,3)) display.brewer.pal(n=20, name="Set1") display.brewer.pal(n=20, name="Set2") display.brewer.pal(n=20, name="Set3") color.gen=c("#F781BF", "#984EA3", "#FF7F00","#66C2A5", "#FC8D62", "#999999","#E78AC3", "#E41A1C", "#8DA0CB", "#A6D854","#FFFF33","#A65628", "#4DAF4A","#377EB8","#BEBADA","#FFD92F", "#E5C494", "#B3B3B3", "#66C2A5", "#FDB462" ) color.gen.trans=c("#F781BF35", "#984EA335", "#FF7F0035","#66C2A535", "#FC8D6235", "#99999935","#E78AC335", "#E41A1C35","#8DA0CB35", "#A6D85435","#FFFF3335", "#A6562835", "#4DAF4A35","#377EB835","#BEBADA35","#FFD92F35", "#E5C49435", "#B3B3B335", "#66C2A535", "#FDB46235" ) genus=c( "Amsinckia","Arabidopsis","Capsella","Clarkia","Collinsia","Dalechampia","Downingia","Erodium","Lasthenia","Leptosiphon","Limnanthes","Medicago","Mimulus","Oenothera sect. anogra","Oenothera sect. oenothera","Polemonium","Primula","Saltugilia","Schiedea","Schizanthus") ########################################## ############ step 6. does range size vary by mating system? ######## ########################################## ### NOTE: first need to save "r squared fuction" at bottom of script to memory (lines 407-end) mixed(log(range.size.res0.05) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(log(range.size.res0.05) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) o=summary(mod)$coefficients[1,1] #grab slope s= summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #grab slope plus intercept exp(o) # back-transformed predicted value outcrossers exp(s) # back-transformed predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus mixed(log(range.size.res0.1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(log(range.size.res0.1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) o=summary(mod)$coefficients[1,1] #grab slope s= summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #grab slope plus intercept exp(o) # back-transformed predicted value outcrossers exp(s) # back-transformed predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus mixed(log(range.size.res0.5) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(log(range.size.res0.5) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) o=summary(mod)$coefficients[1,1] #grab slope s= summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #grab slope plus intercept exp(o) # back-transformed predicted value outcrossers exp(s) # back-transformed predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus mixed(log(range.size.res1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(log(range.size.res1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) o=summary(mod)$coefficients[1,1] #grab slope s= summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #grab slope plus intercept exp(o) # back-transformed predicted value outcrossers exp(s) # back-transformed predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus ########################################## ############ Step 7. Figure 1 - Range size by mating system ######## ########################################## par(mfrow=c(1,1)) # One figure in row 1 and two figures in row 2 # row 1 is 1/3 the height of row 2 # column 2 is 1/4 the width of the column 1 attach(mtcars) layout(matrix(c(1,1,2,2), 2, 2, byrow = TRUE),    widths=c(1,1), heights=c(2,1)) library(arm) tmp3$mate.sys = factor(tmp3$mate.sys,c("selfer","outcrosser")) #slopes allowed to vary across genera (cant include pairs nested within genus, not enough d.f.) mod=lmer(log(range.size.res0.1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) slopes.gen=coef(mod)$genus colnames(slopes.gen)=c("intercept","slope") tmp4=cbind(tmp3,fitted(mod)) par(cex.lab=1.5) # Add extra space to right of plot area; change clipping to figure par(mar=c(3.1, 5.1, 1.1, 14.1), xpd=TRUE) #par(cex.axis=1.5) tmp4$mate.sys = factor(tmp4$mate.sys,c("selfer","outcrosser")) wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("range size km"^"2")), main="",outcol="white",ylim=c(4,12),names=c("",""), boxwex=0.5,border = c("black","black"),yaxt='n', cex.lab=1.7) mtext("selfer", side=1, line=1, at=1,cex=1.6) mtext("outcrosser", side=1, line=1, at=2, cex=1.6) axis(side=2, at=c(5.3,6.9,8.517,10.086,11.84), labels=c("200","1000","5000","24,000","140,000")) for (i in 1:nrow(slopes.gen)) { gen = rownames(slopes.gen)[i] color=color.gen[i] segments(1,slopes.gen$intercept+slopes.gen$slope,2,slopes.gen$intercept, col=color.gen, lwd=0.7, lty=1) } legend("topright", inset=c(-0.6,-0.05), legend=genus, pch=1 , lwd=1,col=color.gen,title="", bty="n", text.font=3) #dev.off() wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("")), main="",outcol="white",ylim=c(4,12),names=c("",""), boxwex=0.5,border = c("black","black"),yaxt='n', add=TRUE, cex.lab=1.7) #### now add differnce barcharts #### #slopes allowed to vary across genera (cant include pairs nested within genus, not enough d.f.) mixed(log(range.size.res0.1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) mod=lmer(log(range.size.res0.1) ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) gen.mean=coef(mod) gen.sd=se.coef(mod) diff=(gen.mean$genus[,1] - (gen.mean$genus[,1]+gen.mean$genus[,2]))*-1 name=rownames(gen.mean$genus) se=gen.sd$genus[,2] diff se par(mar=c(1.1, 5.1, 1.1, 1.1), xpd=TRUE) mp=barplot(diff, axes=FALSE, axisnames=TRUE, col=color.gen,xlab="", ylab="Difference: selfer - outcrosser", ylim=c(-1,2.5),cex.lab=0.9) #box() axis(2,at=seq(-1,2.5,by=1.0)) #abline(h=0) segments(mp,diff-se,mp,diff+se,lwd=1) segments(-0.8,0,25,0) ########################################## ############ Step 8. Does latitude vary by mating system? ######## ########################################## tmp3$lat.mid = ((tmp3$lat.max - tmp3$lat.min) / 2) + tmp3$lat.min mixed(lat.min ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(lat.min ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) summary(mod)$coefficients[1,1] #predicted value outcrossers summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus # mixed(lat.mid ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(lat.mid ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) summary(mod)$coefficients[1,1] #predicted value outcrossers summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus mixed(lat.max ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3,method="LRT") mod=lmer(lat.max ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) summary(mod)$coefficients[1,1] #predicted value outcrossers summary(mod)$coefficients[1,1] + summary(mod)$coefficients[2,1] #predicted value selfers t=list(mod) rsquared.glmm(list(mod)) #proportion of variance explained by mating system and random effects (marginal,conditional) coef(mod)$genus ########################################## ############ Step 9. Figure 2 -- Latitude by mating system ######## ########################################## ####### max ######## mod=lmer(lat.max ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) slopes.gen=coef(mod)$genus colnames(slopes.gen)=c("intercept","slope") tmp4=cbind(tmp3,fitted(mod)) m <- rbind(c(1, 2,3), c(4, 5,6)) layout(m,widths=c(1,1,1), heights=c(2,1) ) # Add extra space to right of plot area; change clipping to figure par(mar=c(3.1, 5.1, 1.1, 1.1), xpd=TRUE) tmp4$mate.sys = factor(tmp4$mate.sys,c("selfer","outcrosser")) wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("maximum latitude")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), ylim=c(10,70),cex.lab=1.6) mtext("selfer", side=1, line=1, at=1,cex=1) mtext("outcrosser", side=1, line=1, at=2, cex=1) mtext("*", side=3, line=-2.5, at=1.5, adj=0.5,cex=1.5) for (i in 1:nrow(slopes.gen)) { gen = rownames(slopes.gen)[i] color=color.gen[i] segments(1,slopes.gen$intercept+slopes.gen$slope,2,slopes.gen$intercept, col=color.gen, lwd=0.7, lty=1) } wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), add=TRUE, ylim=c(10,70),cex.lab=1.6) #### now add differnce barcharts #### gen.mean=coef(mod) gen.sd=se.coef(mod) diff.max=(gen.mean$genus[,1] - (gen.mean$genus[,1]+gen.mean$genus[,2]))*-1 name=rownames(gen.mean$genus) se.max=gen.sd$genus[,2] ####### midpoint ######## mod=lmer(lat.mid ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) slopes.gen=coef(mod)$genus colnames(slopes.gen)=c("intercept","slope") tmp4=cbind(tmp3,fitted(mod)) # Add extra space to right of plot area; change clipping to figure par(mar=c(3.1, 5.1, 1.1, 1.1), xpd=TRUE) tmp4$mate.sys = factor(tmp4$mate.sys,c("selfer","outcrosser")) wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("midpoint latitude")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), ylim=c(10,70),cex.lab=1.6) mtext("selfer", side=1, line=1, at=1,cex=1) mtext("outcrosser", side=1, line=1, at=2, cex=1) mtext("NS", side=3, line=-2.5, at=1.5, adj=0.5,cex=1) for (i in 1:nrow(slopes.gen)) { gen = rownames(slopes.gen)[i] color=color.gen[i] segments(1,slopes.gen$intercept+slopes.gen$slope,2,slopes.gen$intercept, col=color.gen, lwd=0.7, lty=1) } wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), add=TRUE, ylim=c(10,70),cex.lab=1.6) #### now add differnce barcharts #### gen.mean=coef(mod) gen.sd=se.coef(mod) diff.mid=(gen.mean$genus[,1] - (gen.mean$genus[,1]+gen.mean$genus[,2]))*-1 name=rownames(gen.mean$genus) se.mid=gen.sd$genus[,2] ####### minimum ######## mod=lmer(lat.min ~ mate.sys + (1+mate.sys|genus) + (1|ba), weights=pp,data=tmp3) slopes.gen=coef(mod)$genus colnames(slopes.gen)=c("intercept","slope") tmp4=cbind(tmp3,fitted(mod)) par(cex.lab=1.5) wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("minimum latitude")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), ylim=c(10,70),cex.lab=1.6) mtext("selfer", side=1, line=1, at=1,cex=1) mtext("outcrosser", side=1, line=1, at=2, cex=1) mtext("NS", side=3, line=-2.5, at=1.5, adj=0.5,cex=1) for (i in 1:nrow(slopes.gen)) { gen = rownames(slopes.gen)[i] color=color.gen[i] segments(1,slopes.gen$intercept+slopes.gen$slope,2,slopes.gen$intercept, col=color.gen, lwd=0.7, lty=1) } wtd.boxplot(tmp4$'fitted(mod)'~tmp4$mate.sys,weights=tmp4$pp,ylab=expression(paste("")), main="",outcol="white", names=c("",""), boxwex=0.5,border = c("black","black"), add=TRUE, ylim=c(10,70),cex.lab=1.6) #### now add differnce barcharts #### gen.mean=coef(mod) gen.sd=se.coef(mod) diff.min=(gen.mean$genus[,1] - (gen.mean$genus[,1]+gen.mean$genus[,2]))*-1 name=rownames(gen.mean$genus) se.min=gen.sd$genus[,2] ####barcharts by clade mp=barplot(diff.max, axes=TRUE, axisnames=TRUE, col=color.gen,xlab="", ylab="Difference: selfer - outcrosser",cex.lab=0.9, ylim=c(-7,7)) #box() #abline(h=0) segments(mp,diff.max-se.max,mp,diff.max+se.max,lwd=1) segments(-0.8,0,25,0) mp=barplot(diff.mid, axes=TRUE, axisnames=TRUE, col=color.gen,xlab="", ylab="",cex.lab=0.9, ylim=c(-7,7)) #box() #abline(h=0) segments(mp,diff.mid-se.mid,mp,diff.mid+se.mid,lwd=1) segments(-0.8,0,25,0) mp=barplot(diff.min, axes=TRUE, axisnames=TRUE, col=color.gen,xlab="", ylab="",cex.lab=0.9, ylim=c(-7,7)) #box() #abline(h=0) segments(mp,diff.min-se.min,mp,diff.min+se.min,lwd=1) segments(-0.8,0,25,0) ########################################## ############ Step 10. Does range size vary by age, mating system, and their interaction? ######## ########################################## mixed(log(range.size.res0.1) ~ log(dist)*mate.sys + (1|ba) + (1+mate.sys|genus) + (1+log(dist)|genus), weights=pp,data=tmp3,method="LRT") mixed(log(range.size.res0.5) ~ log(dist)*mate.sys + (1|ba) + (1+mate.sys|genus) + (1+log(dist)|genus), weights=pp,data=tmp3,method="LRT") #MODEL FAILS TO CONVERGE! mixed(log(range.size.res1) ~ log(dist)*mate.sys + (1|ba) + (1+mate.sys|genus) + (1+log(dist)|genus), weights=pp,data=tmp3,method="LRT") mixed(log(range.size.res0.05) ~ log(dist)*mate.sys + (1|ba) + (1+mate.sys|genus) + (1+log(dist)|genus), weights=pp,data=tmp3,method="LRT") ########################################## ############ Step 11. Figure 3 -- range size by age, mating system, and their interaction ######## ########################################## par(mfrow=c(1,1), mar=c(5.1,5.1,1,1)) mod=lmer(log(range.size.res0.1) ~ log(dist)*mate.sys + (1|ba) + (1+mate.sys|genus) + (1+log(dist)|genus),, weights=pp,data=tmp3) coef(mod) tmp4=cbind(tmp3,fitted(mod)) selfs = subset(tmp4, mate.sys=="selfer") outs = subset(tmp4, mate.sys=="outcrosser") plot(selfs$'fitted(mod)'~ log(selfs$dist), col="white", ylab=expression(paste("range size km"^"2")), main="",xlab="divergence time (mya)",cex= selfs$fig.weight, bty="n", ylim=c(4,12),yaxt='n',xaxt='n', cex.lab=1.5, cex.main=1.5) axis(side=2, at=c(4,5.99,8.006,9.9988,11.98), labels=c("55","400","3000","22,000","160,000"), cex.axis=0.8) axis(side=1, at=c(-2,-0.9162907,0,1.098612,2.079442,3,4), labels=c("0.14","0.4","1","3","8","20", "55"), cex.axis=0.8) for (i in 1:length(unique(outs$genus))) { gen=sort(unique(outs$genus))[i] color=color.gen[i] tmp.gen = subset(outs, genus==gen) points(tmp.gen$'fitted(mod)' ~ log(dist), data=tmp.gen, col="darkgray", cex= tmp.gen$fig.weight*2, lwd=1.5, pch=1) } for (i in 1:length(unique(selfs$genus))) { gen=sort(unique(selfs$genus))[i] color=color.gen[i] tmp.gen = subset(selfs, genus==gen) points(tmp.gen$'fitted(mod)' ~ log(dist), data=tmp.gen, col="black", cex= tmp.gen$fig.weight*2, lwd=1.5, pch=1) } legend("topright",legend=c("selfer","outcrosser"),pch=1,col=c("black","darkgray"),cex=0.9,lwd=1.3) abline(lm(selfs$'fitted(mod)' ~log(selfs$dist),weights=selfs$pp),col="black") abline(lm(outs$'fitted(mod)' ~log(outs$dist),weights=outs$pp),col="darkgray") ########################################## ############ Step 12. Accounting for correlated traits -- remove pairs with variable ploidy and life history ######## ########################################## sis <-read.csv ("sister.summary.table.csv") sis.sub=subset(sis,mate.diff=="s_o") sis.sub$mate.diff <- factor(sis.sub$mate.diff) sis.sub$ploidy.diff = sis.sub$ploidy.spA - sis.sub$ploidy.spB sis.sub$life.history.diff = sis.sub$life.history.spA - sis.sub$life.history.spB #subset date to just sister pairs with same ploidy and same life.history sis.sub.ploidy=subset(sis.sub, ploidy.diff==0) sis.sub.life=subset(sis.sub, life.history.diff==0) # now run step 4 - 11 on the ploidy and life history data subsets seperately ########################################## ############ r squared functions ######## ########################################## rsquared.glmm <- function(modlist) { # Iterate over each model in the list do.call(rbind, lapply(modlist, r.squared)) }   #' R-squared and pseudo-rsquared for (generalized) linear (mixed) models #' #' This generic function calculates the r squared and pseudo r-squared for #' a variety of(generalized) linear (mixed) model fits. #' Currently implemented for \code{\link{lm}}, \code{\link{lmerTest::merMod}}, #' and \code{\link{nlme::lme}} objects. #' Implementing methods usually call \code{\link{.rsquared.glmm}} #' #' @param mdl a fitted (generalized) linear (mixed) model object #' @return Implementing methods usually return a dataframe with "Class", #' "Family", "Marginal", "Conditional", and "AIC" columns r.squared <- function(mdl){ UseMethod("r.squared") }   #' Marginal r-squared for lm objects #' #' This method uses r.squared from \code{\link{summary}} as the marginal. #' Contrary to other \code{\link{r.squared}} methods, #' this one doesn't call \code{\link{.rsquared.glmm}} #' #' @param mdl an lm object (usually fit using \code{\link{lm}}, #' @return a dataframe with with "Class" = "lm", "Family" = "gaussian", #' "Marginal" = unadjusted r-squared, "Conditional" = NA, and "AIC" columns r.squared.lm <- function(mdl){ data.frame(Class=class(mdl), Family="gaussian", Link="identity", Marginal=summary(mdl)$r.squared, Conditional=NA, AIC=AIC(mdl)) }   #' Marginal and conditional r-squared for merMod objects #' #' This method extracts the variance for fixed and random effects, residuals, #' and the fixed effects for the null model (in the case of Poisson family), #' and calls \code{\link{.rsquared.glmm}} #' #' @param mdl an merMod model (usually fit using \code{\link{lme4::lmer}}, #' \code{\link{lme4::glmer}}, \code{\link{lmerTest::lmer}}, #' \code{\link{blme::blmer}}, \code{\link{blme::bglmer}}, etc) r.squared.merMod <- function(mdl){ # Get variance of fixed effects by multiplying coefficients by design matrix VarF <- var(as.vector(lme4::fixef(mdl) %*% t(mdl@pp$X))) # Get variance of random effects by extracting variance components # Omit random effects at the observation level, variance is factored in later VarRand <- sum( sapply( VarCorr(mdl)[!sapply(unique(unlist(strsplit(names(ranef(mdl)),":|/"))), function(l) length(unique(mdl@frame[,l])) == nrow(mdl@frame))], function(Sigma) { X <- model.matrix(mdl) Z <- X[,rownames(Sigma)] sum(diag(Z %*% Sigma %*% t(Z)))/nrow(X) } ) ) # Get the dispersion variance VarDisp <- unlist(VarCorr(mdl)[sapply(unique(unlist(strsplit(names(ranef(mdl)),":|/"))), function(l) length(unique(mdl@frame[,l])) == nrow(mdl@frame))]) if(is.null(VarDisp)) VarDisp = 0 else VarDisp = VarDisp if(inherits(mdl, "lmerMod")){ # Get residual variance VarResid <- attr(lme4::VarCorr(mdl), "sc")^2 # Get ML model AIC mdl.aic <- AIC(update(mdl, REML=F)) # Model family for lmer is gaussian family <- "gaussian" # Model link for lmer is identity link <- "identity" } else if(inherits(mdl, "glmerMod")){ # Get the model summary mdl.summ <- summary(mdl) # Get the model's family, link and AIC family <- mdl.summ$family link <- mdl.summ$link mdl.aic <- AIC(mdl) # Pseudo-r-squared for poisson also requires the fixed effects of the null model if(family=="poisson") { # Get random effects names to generate null model rand.formula <- reformulate(sapply(findbars(formula(mdl)), function(x) paste0("(", deparse(x), ")")), response=".") # Generate null model (intercept and random effects only, no fixed effects) null.mdl <- update(mdl, rand.formula) # Get the fixed effects of the null model null.fixef <- as.numeric(lme4::fixef(null.mdl)) } } # Call the internal function to do the pseudo r-squared calculations .rsquared.glmm(VarF, VarRand, VarResid, VarDisp, family = family, link = link, mdl.aic = mdl.aic, mdl.class = class(mdl), null.fixef = null.fixef) }   #' Marginal and conditional r-squared for lme objects #' #' This method extracts the variance for fixed and random effects, #' as well as residuals, and calls \code{\link{.rsquared.glmm}} #' #' @param mdl an lme model (usually fit using \code{\link{nlme::lme}}) r.squared.lme <- function(mdl){ # Get design matrix of fixed effects from model Fmat <- model.matrix(eval(mdl$call$fixed)[-2], mdl$data) # Get variance of fixed effects by multiplying coefficients by design matrix VarF <- var(as.vector(nlme::fixef(mdl) %*% t(Fmat))) # Get variance of random effects by extracting variance components VarRand <- sum(suppressWarnings(as.numeric(nlme::VarCorr(mdl) [rownames(nlme::VarCorr(mdl)) != "Residual", 1])), na.rm=T) # Get residual variance VarResid <- as.numeric(nlme::VarCorr(mdl)[rownames(nlme::VarCorr(mdl))=="Residual", 1]) # Call the internal function to do the pseudo r-squared calculations .rsquared.glmm(VarF, VarRand, VarResid, family = "gaussian", link = "identity", mdl.aic = AIC(update(mdl, method="ML")), mdl.class = class(mdl)) }   #' Marginal and conditional r-squared for glmm given fixed and random variances #' #' This function is based on Nakagawa and Schielzeth (2013). It returns the marginal #' and conditional r-squared, as well as the AIC for each glmm. #' Users should call the higher-level generic "r.squared", or implement a method for the #' corresponding class to get varF, varRand and the family from the specific object #' #' @param varF Variance of fixed effects #' @param varRand Variance of random effects #' @param varResid Residual variance. Only necessary for "gaussian" family #' @param family family of the glmm (currently works with gaussian, binomial and poisson) #' @param link model link function. Working links are: gaussian: "identity" (default); #' binomial: "logit" (default), "probit"; poisson: "log" (default), "sqrt" #' @param mdl.aic The model's AIC #' @param mdl.class The name of the model's class #' @param null.fixef Numeric vector containing the fixed effects of the null model. #' Only necessary for "poisson" family #' @return A data frame with "Class", "Family", "Marginal", "Conditional", and "AIC" columns .rsquared.glmm <- function(varF, varRand, varResid = NULL, varDisp = NULL, family, link, mdl.aic, mdl.class, null.fixef = NULL){ if(family == "gaussian"){ # Only works with identity link if(link != "identity") family_link.stop(family, link) # 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) } else if(family == "binomial"){ # Get the distribution-specific variance if(link == "logit") varDist <- (pi^2)/3 else if(link == "probit") varDist <- 1 else family_link.stop(family, link) # Calculate marginal R-squared Rm <- varF/(varF+varRand+varDist+varDisp) # Calculate conditional R-squared (fixed effects+random effects/total variance) Rc <- (varF+varRand)/(varF+varRand+varDist+varDisp) } else if(family == "poisson"){ # Get the distribution-specific variance if(link == "log") varDist <- log(1+1/exp(null.fixef)) else if(link == "sqrt") varDist <- 0.25 else family_link.stop(family, link) # Calculate marginal R-squared Rm <- varF/(varF+varRand+varDist+varDisp) # Calculate conditional R-squared (fixed effects+random effects/total variance) Rc <- (varF+varRand)/(varF+varRand+varDist+varDisp) } else family_link.stop(family, link) # Bind R^2s into a matrix and return with AIC values data.frame(Class=mdl.class, Family = family, Link = link, Marginal=Rm, Conditional=Rc, AIC=mdl.aic) }   #' stop execution if unable to calculate variance for a given family and link family_link.stop <- function(family, link){ stop(paste("Don't know how to calculate variance for", family, "family and", link, "link.")) }