Code to Re-Analyze the results of Maestre et al 2012 Spurred on by Bob O'Hara, Jeremy Fox, and Ted Hart's ongoing discussions

Maestre_reanalysis.R

########################################################################################################
# Re-Analysis of Maestre et al 2012
# examining relationship between biodiversity and ecosystem multifunctionality
# Paper and data included as a supplement can be found at  http://dx.doi.org/10.1126/science.1215442
########################################################################################################

######################
####### First, let's look at standardized regression coefficients
####### Both for the best fit model, but also for the model averaged coefficients
######################

library(QuantPsyc)
library(AICcmodavg)

mdata<-read.csv("./Maestre_2012_1215442databaseS1.csv")
pairs(mdata[,4:17]) #plotting for fun and profit

#fit the series of models based on Table 2
mod1<-lm(MUL ~ SR+SLO+SAC+PCA_C1+PCA_C3+PCA_C4+LAT+LONG+ELE, data=mdata)
mod2<-update(mod1, . ~ . -PCA_C3)
mod3<-update(mod2, . ~ . -PCA_C1)
mod4<-update(mod1, . ~ . -PCA_C3 - PCA_C4)
mod5<-update(mod1, . ~ . +PCA_C2)
mod6<-update(mod1, . ~ . -PCA_C4)
mod7<-update(mod5, . ~ . -PCA_C1)
mod8<-update(mod5, . ~ . -PCA_C3)

#AIC table - results slightly different - what is wrong above?
aictab(list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8), 1:8)

#look at unstandardized and standardized coefficients
summary(mod1)
lm.beta(mod1)

#what about the standardized effect from the averaged model?
library(AICcmodavg)
srEffect<-modavg(cand.set =list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8), parm="SR", modnames=1:8)$Mod.avg.beta

srEffect*sd(mdata$SR, na.rm=T)/sd(mdata$MUL, na.rm=T)


#look at all standardized coefficients from the averaged model
sapply(names(mdata)[4:13], function(x) {
  Effect<-modavg(cand.set =list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8), parm=x, modnames=1:8)$Mod.avg.beta
  Effect*sd(mdata[[x]], na.rm=T)/sd(mdata$MUL, na.rm=T)
})


####################
### Part 2 - Is the diversity effect just reflecting
### diversity being caused by the same predictors as
### multifunction?  An SEM comparison
####################

library(lavaan)

#Model 1: Multifunction driven by species richness and other predictors
model1<-'
  SR ~ SLO+SAC+PCA_C1+PCA_C3+PCA_C4+LAT+LONG+ELE
  MUL ~ SR+SLO+SAC+PCA_C1+PCA_C3+PCA_C4+LAT+LONG+ELE
'
fit1<-sem(model1, data=mdata)

#Model 2: Richness & Multifunction driven by the same causes
# Note: We're assuming that this is a sufficient set of drivers
# as otherwise if SR and MUL covary, this is an equivalent model to model 1
model2<-'
  SR ~ SLO+SAC+PCA_C1+PCA_C3+PCA_C4+LAT+LONG+ELE
  MUL ~ SLO+SAC+PCA_C1+PCA_C3+PCA_C4+LAT+LONG+ELE
  SR ~~ 0*MUL #assumption that there is no correlation beyond the set of common drivers
'


fit2<-sem(model2, data=mdata)

#compare the two models with a Likelihood Ratio Test
anova(fit1, fit2)


summary(fit1, standardized=T, rsquare=T)