#These code for the statistical freeware package R were used in the data analyses in ‘Bradford MA, Veen, G.F.C., Bonis, A., Bradford, E.M., Classen, A.T., Cornelissen, J.H.C., Crowther, T.W., De Long, J.R., Freschet, G.T., Kardol, P., Manrubia-Freixa, M., Maynard, D.S., Newman, G.S., van Logtestijn, R.S.P., Viketoft, M., Wardle, D.A., Wieder, W.R., Wood, S.A., van der Putten, W.H. Testing the hierarchical model of litter decomposition. Nature Ecology and Evolution, in press.’


#Data and code, along with a ReadMe file (BradfordVeenLitDecompREADME.txt), were submitted to the Dryad data repository August 24th 2017 following the papers acceptance.

#The persons associated with this code have dedicated the work to the public domain by waiving all of their rights to the code worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. You can copy, modify, distribute and use the code, even for commercial purposes, all without asking permission.

#Code prepared by Mark A Bradford, Yale University. mark.bradford@yale.edu

#Read in the individual observations and potential explanatory variables. Data are organized to look at effects on litter decomposition (Holcus and Festuca)

read.csv("/Users/mab25/Desktop/BradfordVeenLitDecompDATA.csv",header=T,na.strings=c("NA"))
Grad2<-read.csv("/Users/mab25/Desktop/BradfordVeenLitDecompDATA.csv",header=T,na.strings=c("NA"))

options(contrasts=c("contr.treatment","contr.poly"))
#install.packages("package name")
library(lme4)
library(MuMIn)
library(sampling)
library(LMERConvenienceFunctions)
library(car)
library(MASS)
library(usdm)
library(arm)
#library(LanguageR)

#Establish factors to account for potential spatial autocorrelation using random effects
Grad2$SITE2<-factor(Grad2$Site2)
is.factor(Grad2$SITE2)
Grad2$TRANSECT<-factor(Grad2$Transect)
is.factor(Grad2$TRANSECT)
Grad2$QUADRAT<-factor(Grad2$Quadrat)
is.factor(Grad2$QUADRAT)

#Check Lit variables are discrete factors
is.factor(Grad2$Lit)
is.factor(Grad2$LitNum)


#calculate a mean of the three spot measures of soil temperature and moisture across the 3 month field incubations
Grad2$meanT<-((Grad2$Temp1+Grad2$Temp2+Grad2$Temp3)/3) # using Grad2$ associates this new variable with the read in dataset Grad2.
Grad2$meanM<-((Grad2$Moist1+Grad2$Moist2+Grad2$Moist3)/3)

#used litter moisture in the final models given site differences in soil texture

#create soil N availability index (initial extractable inorganic N)
Grad2$Navail<-(Grad2$NOxinit+Grad2$NH4init)

#r2 function for evaluating models:

r2.mixed<-function(mF){
  mFX<-model.matrix(mF)
  VarF <- var(as.vector(fixef(mF) %*% t(mFX)))
  VarR<-sum(as.numeric(VarCorr(mF))) 
  VarResid<-attr(VarCorr(mF), "sc")^2
  fR2<-VarF/(VarF + VarR + VarResid)
  rfR2<-(VarF + VarR)/(VarF + VarR + VarResid)
  list(fR2=fR2,rfR2=rfR2)
}


#vif function for evaluating models:

vif.mer <- function (fit) {
  ## adapted from rms::vif
  
  v <- vcov(fit)
  nam <- names(fixef(fit))
  
  ## exclude intercepts
  ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
  if (ns > 0) {
    v <- v[-(1:ns), -(1:ns), drop = FALSE]
    nam <- nam[-(1:ns)]
  }
  
  d <- diag(v)^0.5
  v <- diag(solve(v/(d %o% d)))
  names(v) <- nam
  v
}

kappa.mer <- function (fit,
                       scale = TRUE, center = FALSE,
                       add.intercept = TRUE,
                       exact = FALSE) {
  X <- fit@X
  nam <- names(fixef(fit))
  
  ## exclude intercepts
  nrp <- sum(1 * (nam == "(Intercept)"))
  if (nrp > 0) {
    X <- X[, -(1:nrp), drop = FALSE]
    nam <- nam[-(1:nrp)]
  }
  
  if (add.intercept) {
    X <- cbind(rep(1), scale(X, scale = scale, center = center))
    kappa(X, exact = exact)
  } else {
    kappa(scale(X, scale = scale, center = scale), exact = exact)
  }
}

colldiag.mer <- function (fit,
                          scale = TRUE, center = FALSE,
                          add.intercept = TRUE) {
  ## adapted from perturb::colldiag, method in Belsley, Kuh, and
  ## Welsch (1980). look for a high condition index (> 30) with
  ## more than one high variance propotion. see ?colldiag for more
  ## tips.
  result <- NULL
  if (center)
    add.intercept <- FALSE
  if (is.matrix(fit) || is.data.frame(fit)) {
    X <- as.matrix(fit)
    nms <- colnames(fit)
  }
  else if (class(fit) == "mer") {
    nms <- names(fixef(fit))
    X <- fit@X
    if (any(grepl("(Intercept)", nms))) {
      add.intercept <- FALSE
    }
  }
  X <- X[!is.na(apply(X, 1, all)), ]
  
  if (add.intercept) {
    X <- cbind(1, X)
    colnames(X)[1] <- "(Intercept)"
  }
  X <- scale(X, scale = scale, center = center)
  
  svdX <- svd(X)
  svdX$d
  condindx <- max(svdX$d)/svdX$d
  dim(condindx) <- c(length(condindx), 1)
  
  Phi = svdX$v %*% diag(1/svdX$d)
  Phi <- t(Phi^2)
  pi <- prop.table(Phi, 2)
  colnames(condindx) <- "cond.index"
  if (!is.null(nms)) {
    rownames(condindx) <- nms
    colnames(pi) <- nms
    rownames(pi) <- nms
  } else {
    rownames(condindx) <- 1:length(condindx)
    colnames(pi) <- 1:ncol(pi)
    rownames(pi) <- 1:nrow(pi)
  }
  
  result <- data.frame(cbind(condindx, pi))
  zapsmall(result)
}

maxcorr.mer <- function (fit,
                         exclude.intercept = TRUE) {
  so <- summary(fit)
  corF <- so@vcov@factors$correlation
  nam <- names(fixef(fit))
  
  ## exclude intercepts
  ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
  if (ns > 0 & exclude.intercept) {
    corF <- corF[-(1:ns), -(1:ns), drop = FALSE]
    nam <- nam[-(1:ns)]
  }
  corF[!lower.tri(corF)] <- 0
  maxCor <- max(corF)
  minCor <- min(corF)
  if (abs(maxCor) > abs(minCor)) {
    zapsmall(maxCor)
  } else {
    zapsmall(minCor)
  }
}


#Subset Holcus and Festuca data
dataHol <- Grad2[Grad2$Lit=="Hol",]
dataFes <- Grad2[Grad2$Lit=="Fes",]

#Plotting quadrat-level explanatory variables for the points by site:
#Using the subsetted Holcus data to avoid duplicating observations per quadrat (as they’re the same for each litter pair)

plot(dataHol$meanT~dataHol$Site2)
plot(dataHol$mLitMoist~dataHol$Site2)
plot(dataHol$SIR~dataHol$Site2)
plot(dataHol$TpH~dataHol$Site2)
plot(dataHol$Navail~dataHol$Site2)

#Getting the ranges in the explanatory variables to generate regression relationships and evaluate them

range(dataHol$meanT,na.rm = TRUE)
range(dataHol$mLitMoist,na.rm = TRUE)
range(dataHol$SIR,na.rm = TRUE)
range(dataHol$TpH,na.rm = TRUE)
range(dataHol$Navail,na.rm = TRUE)

————————————————————————————————————
————————————————————————————————————

#ANALYZING FULL DATASET FOR LITTER C LOSS

#plotting the response variable distributions to look at assumptions of normality:
hist(Grad2$SLitCloss)
qqnorm(Grad2$SLitCloss)
qqline(Grad2$SLitCloss)

#(1|SITE2:TRANSECT:QUADRAT) NESTS THE RANDOM EFFECTS

LitCmodel1<-lmer(SLitCloss~LitNum+meanT+mLitMoist+Navail+TpH+SIR+LitNum:meanT+LitNum:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel1)
summary(standardize(LitCmodel1))
pamer.fnc(LitCmodel1)
plot(LitCmodel1)
r2.mixed(LitCmodel1)
sqrt(vif.mer(standardize(LitCmodel1)))
#want values for each fixed term that are <2

#SCALING THE RESPONSE VARIABLE SO CAN COMPARE TO WITHIN-SITE MODELS

LitCmodel2<-lmer(scale(SLitCloss)~LitNum+meanT+mLitMoist+Navail+TpH+SIR+LitNum:meanT+LitNum:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel2)
summary(standardize(LitCmodel2))
pamer.fnc(LitCmodel2)
plot(LitCmodel2)
r2.mixed(LitCmodel2)
sqrt(vif.mer(standardize(LitCmodel2)))

#Checking for possibility of nonlinear responses to temperature and moisture (strong evidence to fit a second order term for temperature)

LitCmodel3<-lmer(scale(SLitCloss)~LitNum+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitNum:meanT+LitNum:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel3)
summary(standardize(LitCmodel3))
pamer.fnc(LitCmodel3)
plot(LitCmodel3)
r2.mixed(LitCmodel3)
sqrt(vif.mer(standardize(LitCmodel3)))


#Testing with Litter Quality represented as initial %N, instead of a binary variable

LitCmodel4<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel4)
summary(standardize(LitCmodel4))
pamer.fnc(LitCmodel4)
plot(LitCmodel4)
r2.mixed(LitCmodel4)
sqrt(vif.mer(standardize(LitCmodel4)))

#Checking that using “scale” to center and divide by 1 SD the response variable does exactly that

LitCmodel5<-lmer(scale(SLitCloss)~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel5)
summary(standardize(LitCmodel5))
pamer.fnc(LitCmodel5)
plot(LitCmodel5)
r2.mixed(LitCmodel5)
sqrt(vif.mer(standardize(LitCmodel5)))

#scaling the response variable

Grad2$CLossSTD <- ((Grad2$SLitCloss)-(mean(Grad2$SLitCloss,na.rm=T)))/(sd(Grad2$SLitCloss,na.rm=T))
mean(Grad2$CLossSTD,na.rm=T)
sd(Grad2$CLossSTD,na.rm=T)

LitCmodel5b<-lmer(CLossSTD~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2)
summary(LitCmodel5b)
summary(standardize(LitCmodel5b))
pamer.fnc(LitCmodel5b)
plot(LitCmodel5b)
r2.mixed(LitCmodel5b)
sqrt(vif.mer(standardize(LitCmodel5b)))

#checking for highly influential/outlier observations: using a linear model given improved model checking

LitCmodel4a<-lm(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail,data=Grad2)
summary(LitCmodel4a)
plot(LitCmodel4a)

cutoff <- 4/((nrow(mtcars)-length(LitCmodel4a$coefficients)-2)) 
plot(LitCmodel4a, which=4, cook.levels=cutoff)

#observation 157 is a substantial outlier. no other variables with enough influence to remove (157 removed on outlier basis in terms of how it affected model fits)

#“Standardize” does not work with the lmer when omitting observations. Therefore generate standardized values for each of the factors per Gelman 2008. Mean should = 0 and sd = 0.5

Grad2$LitNSTD <- ((Grad2$LitinitN)-(mean(Grad2$LitinitN,na.rm=T)))/(2*sd(Grad2$LitinitN,na.rm=T))
mean(Grad2$LitNSTD)
sd(Grad2$LitNSTD)

Grad2$meanTSTD <- ((Grad2$meanT)-(mean(Grad2$meanT,na.rm=T)))/(2*sd(Grad2$meanT,na.rm=T))
mean(Grad2$meanTSTD,na.rm=T)
sd(Grad2$meanTSTD,na.rm=T)

Grad2$moistSTD <- ((Grad2$mLitMoist)-(mean(Grad2$mLitMoist,na.rm=T)))/(2*sd(Grad2$mLitMoist,na.rm=T))
mean(Grad2$moistSTD,na.rm=T)
sd(Grad2$moistSTD,na.rm=T)

Grad2$NavSTD <- ((Grad2$Navail)-(mean(Grad2$Navail,na.rm=T)))/(2*sd(Grad2$Navail,na.rm=T))
mean(Grad2$NavSTD,na.rm=T)
sd(Grad2$NavSTD,na.rm=T)

Grad2$pHSTD <- ((Grad2$TpH)-(mean(Grad2$TpH,na.rm=T)))/(2*sd(Grad2$TpH,na.rm=T))
mean(Grad2$pHSTD,na.rm=T)
sd(Grad2$pHSTD,na.rm=T)

Grad2$sirSTD <- ((Grad2$SIR)-(mean(Grad2$SIR,na.rm=T)))/(2*sd(Grad2$SIR,na.rm=T))
mean(Grad2$sirSTD,na.rm=T)
sd(Grad2$sirSTD,na.rm=T) 

#Check that get same output whether use “standardize” or calculate the standard values directly

LitCmodel5c<-lm(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail,data=Grad2)
summary(standardize(LitCmodel5c))

LitCmodel5d<-lm(SLitCloss~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+NavSTD+pHSTD+sirSTD+LitNSTD:meanTSTD+LitNSTD:NavSTD,data=Grad2)
summary(LitCmodel5d)

#omitting observation 157

LitCmodel4c<-lm(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail,data=Grad2[-c(157),])
summary(LitCmodel4c)
plot(LitCmodel4c)

cutoff <- 4/((nrow(mtcars)-length(LitCmodel4c$coefficients)-2)) 
plot(LitCmodel4c, which=4, cook.levels=cutoff)

#omitting in the lmer framework

LitCmodel6<-lmer(SLitCloss~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+NavSTD+pHSTD+sirSTD+LitNSTD:meanTSTD+LitNSTD:NavSTD+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel6)
pamer.fnc(LitCmodel6)
plot(LitCmodel6)
r2.mixed(LitCmodel6)
sqrt(vif.mer(LitCmodel6))

#standardizing the response variable (mean of 0, stdev of 1 for the response variable)

Grad2$CLossSTD <- ((Grad2$SLitCloss)-(mean(Grad2$SLitCloss,na.rm=T)))/(sd(Grad2$SLitCloss,na.rm=T))
mean(Grad2$CLossSTD,na.rm=T)
sd(Grad2$CLossSTD,na.rm=T)

LitCmodel7<-lmer(CLossSTD~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+NavSTD+pHSTD+sirSTD+LitNSTD:meanTSTD+LitNSTD:NavSTD+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel7)
pamer.fnc(LitCmodel7)
plot(LitCmodel7)
r2.mixed(LitCmodel7)
sqrt(vif.mer(LitCmodel7))

#LitCmodel8 is the same as model7 but with unstandardized coefficients (note that the vifs for terms associated with litter quality, N availability and temperature are elevated (much greater than 2) when run with non-standardized variables)

LitCmodel8<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+TpH+SIR+LitinitN:meanT+LitinitN:Navail+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel8)
pamer.fnc(LitCmodel8)
plot(LitCmodel8)
r2.mixed(LitCmodel8)
sqrt(vif.mer(LitCmodel8))

#LitCmodel9 is the same as model7 but with the local variables (i.e. pH, SIR and soil N availability) dropped. VIFs all >2 except for moisture

LitCmodel9<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+LitinitN:meanT+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel9)
pamer.fnc(LitCmodel9)
plot(LitCmodel9)
r2.mixed(LitCmodel9)
sqrt(vif.mer(LitCmodel9))

#LitCmodel10 is the same as model9 but with site mean values for the litter and microclimate factors. VIFs all >2 except for moisture

LitCmodel10<-lmer(SLitCloss~MEANlin+MEANt+I(MEANt^2)+MEANlim+MEANlin:MEANt+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel10)
pamer.fnc(LitCmodel10)
plot(LitCmodel10)
r2.mixed(LitCmodel10)
sqrt(vif.mer(LitCmodel10))

#LitCmodel11 is the same as model9 but with site mean values for the litter and microclimate factors, and the site mean values for the "local" factors. VIFs all >2 except for moisture

LitCmodel11<-lmer(SLitCloss~MEANlin+MEANt+I(MEANt^2)+MEANlim+MEANnava+MEANtph+MEANsir+MEANlin:MEANt+MEANlin:MEANnava+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel11)
pamer.fnc(LitCmodel11)
plot(LitCmodel11)
r2.mixed(LitCmodel11)
sqrt(vif.mer(LitCmodel11))

---------------------------
#Updates to above code

#LitCmodel10b is the same as model15 (see below) but with site mean values for the main factors only in the data set (but maintains the n)

LitCmodel10b<-lmer(SLitCloss~MEANlin+MEANt+I(MEANt^2)+MEANlim+MEANnava+MEANsir+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel10b)
pamer.fnc(LitCmodel10b)
plot(LitCmodel10b)
r2.mixed(LitCmodel10b)
sqrt(vif.mer(LitCmodel10b))

#Dropping the factors not usually measured (SIR and N avail) (REPORTED IN THE MANUSCRIPT AS THE SITE MEAN LQ+CLIMATE MODEL)

LitCmodel10c<-lmer(SLitCloss~MEANlin+MEANt+I(MEANt^2)+MEANlim+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel10c)
pamer.fnc(LitCmodel10c)
plot(LitCmodel10c)
r2.mixed(LitCmodel10c)
sqrt(vif.mer(LitCmodel10c))


---------------------------
#Additional models to assess whether need more interactions included, such as temperature * moisture

#Running first with the standardized coefficients, building from the full model #7 above but dropping pH and focusing on conceptual model of N, litter, SIR and microclimate only. NOTE ALSO THAT THE RESPONSE VARIABLE WAS NOT STANDARDIZED (REPORTED IN THE MANUSCRIPT AS THE MICROSITE WITH INTERACTIONS MODEL: STANDARDIZED COEFFICIENTS)

LitCmodel12<-lmer(SLitCloss~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+NavSTD+sirSTD+LitNSTD:meanTSTD+LitNSTD:moistSTD+LitNSTD:sirSTD+LitNSTD:NavSTD+meanTSTD:moistSTD+meanTSTD:NavSTD+meanTSTD:sirSTD+moistSTD:NavSTD+moistSTD:sirSTD+NavSTD:sirSTD+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel12)
pamer.fnc(LitCmodel12)
plot(LitCmodel12)
r2.mixed(LitCmodel12)
sqrt(vif.mer(LitCmodel12))

#Explains very much the same variation as Model7. Testing whether plotting the unstandardized coefficients markedly changes the outcomes of the plotted regressions vs. Model7 (it doesn't). (REPORTED IN THE MANUSCRIPT AS THE MICROSITE WITH INTERACTIONS MODEL: UNSTANDARDIZED COEFFICIENTS)

LitCmodel13<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+SIR+LitinitN:meanT+LitinitN:mLitMoist+LitinitN:SIR+LitinitN:Navail+meanT:mLitMoist+meanT:Navail+meanT:SIR+mLitMoist:Navail+mLitMoist:SIR+Navail:SIR+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel13)
pamer.fnc(LitCmodel13)
plot(LitCmodel13)
r2.mixed(LitCmodel13)
sqrt(vif.mer(LitCmodel13))

#Dropping all of the interactions to evaluate whether the model gives markedly different R2 and regression outputs (note response variable not standardized)

LitCmodel14<-lmer(SLitCloss~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+NavSTD+sirSTD+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel14)
pamer.fnc(LitCmodel14)
plot(LitCmodel14)
r2.mixed(LitCmodel14)
sqrt(vif.mer(LitCmodel14))

#Unstandardized (REPORTED IN THE MANUSCRIPT AS THE MICROSITE MAIN EFFECT MODEL)

LitCmodel15<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+Navail+SIR+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel15)
pamer.fnc(LitCmodel15)
plot(LitCmodel15)
r2.mixed(LitCmodel15)
sqrt(vif.mer(LitCmodel15))



#Dropping the "local" factors (SIR and N avail) in the main effect only model, keeping variance

LitCmodel16<-lmer(SLitCloss~LitNSTD+meanTSTD+I(meanTSTD^2)+moistSTD+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel16)
pamer.fnc(LitCmodel16)
plot(LitCmodel16)
r2.mixed(LitCmodel16)
sqrt(vif.mer(LitCmodel16))

#Unstandardized (REPORTED IN THE MANUSCRIPT AS THE MICROSITE LQ+CLIMATE MODEL)

LitCmodel17<-lmer(SLitCloss~LitinitN+meanT+I(meanT^2)+mLitMoist+(1|SITE2:TRANSECT:QUADRAT),data=Grad2[-c(157),])
summary(LitCmodel17)
pamer.fnc(LitCmodel17)
plot(LitCmodel17)
r2.mixed(LitCmodel17)
sqrt(vif.mer(LitCmodel17))

—————————————————————
—————————————————————