################################################################### #### Statistics and calculation of maintenance feeding #### for Lang, Ehnes, Brose & Rall 2017: Temperature and consumer type dependencies of energy flows in natural communities #### doi: 10.1111/oik.04419 #### last edited 05.04.2017 #### 2017, the authors ################################################################### #### for the last time performed by Bjoern C. Rall under R-Version: # R version 3.3.3 (2017-03-06) -- "Another Canoe" # Copyright (C) 2017 The R Foundation for Statistical Computing # Platform: x86_64-pc-linux-gnu (64-bit) #### # R Core Team (2017). R: A language and environment for statistical # computing. R Foundation for Statistical Computing, Vienna, Austria. # URL https://www.R-project.org/. ################################################################### #### Packages used and their citations: #### #### nlme (v. 3.1-128) for mixed effects models # Pinheiro J, Bates D, DebRoy S, Sarkar D and R Core Team (2016). _nlme: # Linear and Nonlinear Mixed Effects Models_. R package version 3.1-128, # . #### #### piecewiseSEM (v. 1.2.1) for R-squared analyses: # Lefcheck, Jonathan S. (2015) piecewiseSEM: Piecewise structural # equation modeling in R for ecology, evolution, and systematics. # Methods in Ecology and Evolution. 7(5): 573-579. DOI: # 10.1111/2041-210X.12512 #### #### gtools (v. 3.5.0) for logit und inv.logit # Gregory R. Warnes, Ben Bolker and Thomas Lumley (2015). gtools: # Various R Programming Tools. R package version 3.5.0. # https://CRAN.R-project.org/package=gtools #### #### For all citations above: # These are the automatically generated citations! #### # Find further information in the main Manuscript, including further references on methods used below. #### ####################################################################################################################################### ################################################## Table of contents ################################################################## ####################################################################################################################################### #### 1) Overall settings and commands.....................................................................................l. 53 #### 2) Statistics for assimilation efficiency............................................................................l. 92 #### 3) Statistics for respiration rates..................................................................................l. 264 #### 4) Calculation of maintenance feeding rates..........................................................................l. 423 #### 5) Statistics for assimilation efficiency using the alternative random structure (sensitivity analyses)..............l. 595 #### 6) Statistics for assimilation efficiency using only originally reported fresh weights................................l. 730 ####################################################################################################################################### ################################################## Overall settings and commands ###################################################### ####################################################################################################################################### #### Removes all old lists/objects: rm(list=ls()) #### Loading the libraries needed (you may have to install them beforehand): library("nlme") # here: for mixed effects models library("piecewiseSEM") # here: for estimating the r-squared library("gtools") # here: for logit transformation #### Citations, etc.: #citation() #citation("nlme") #citation("piecewiseSEM") #citation("gtools") #print(sessionInfo()) #### Setting the working directory and load the data: setwd("/homes/br12hony/Dokumente/projects/ZZZ_Other/LangBirgit_2009_Assimilation/2017-04-12") ds.raw = read.csv("Lang_et_al_2017_data.csv") # The data contains assimilation efficiency and respiration data (ds = dataset). str(ds.raw) # Please check if the data was correctly read in. #### Needed transformations: ## The logit transformation is needed to assure linearity and that predictions do not cross 0 and 1, see methods for further information: ds.raw$ae.logit = gtools::logit(ds.raw$assimilation.efficiency,min=0,max=1) ## The ln-transformation to assure linearity and non-negative metabolic rates: ds.raw$log.met = log(ds.raw$metabolism.J.h) ## The ln-transformation to assure linearity and non-negative body masses: ds.raw$log.weight = log(ds.raw$body.size.gram) ## Temperature transformation, see methods and references therein for details: boltz = 0.00008617343 # The Boltzmann constant T0 = 293.15 # The temperature where the estimated constant (intercept in after log transformation) is shifted to, here 20°C <-> 293.15 K. ds.raw$temp.kelvin = 273.15+ds.raw$temperature.degree.C # The temperature in the dataset transformed to Kelvin. ds.raw$temp.kT = (ds.raw$temp.kelvin-T0)/(boltz*ds.raw$temp.kelvin*T0) # The transformation of temperature data to the "Arrhenius temperature". ####################################################################################################################################### ################################################## Statistics for assimilation efficiency ############################################# ####################################################################################################################################### ds.ae = subset(ds.raw, assimilation.efficiency > -999) # Select only data for assimilation efficiency (ds.ae = dataset containing only assimilation efficiency). nrow(ds.ae) # number of data points #### Selection for the best random stucture and heterogeneity - we use REML here, see Methods for details: lme.ae.weights.REML = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.REML = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML" ) lme.ae.weights.REML.b = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.REML.b = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "REML" ) ## Calculating the BIC and the delta-BIC for the 4 REML models with varied random structure and data heterogeneity (Table S3): BIC(lme.ae.weights.REML,lme.ae.REML, lme.ae.weights.REML.b,lme.ae.REML.b) BIC(lme.ae.weights.REML,lme.ae.REML, lme.ae.weights.REML.b,lme.ae.REML.b)[2]-min(BIC(lme.ae.weights.REML,lme.ae.REML, lme.ae.weights.REML.b,lme.ae.REML.b)[2]) ## The model "lme.ae.weights.REML.b" is the best, but less than 2 delta-BIC units apart from the next best model - we choose this random structure with weights for further selection of the fixed effects. ## We also performed a sensitivity analyses on the fixed structure using the random structure of model "lme.ae.weights.REML", see below. #### Selection for the best fixed structure, we use ML for model selection - see methods and references therein for details: # The full model: lme.ae.5 = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding one interaction between explanatory variables: lme.ae.4a = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.4b = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, or one interaction and one explanatory variable: lme.ae.3a = lme(ae.logit ~ log.weight + temp.kT + consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.3b = lme(ae.logit ~ log.weight + consumer.type + consumer.type:log.weight, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.3c = lme(ae.logit ~ temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, and one explanatory variable: lme.ae.2a = lme(ae.logit ~ log.weight + temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.2b = lme(ae.logit ~ log.weight + consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.2c = lme(ae.logit ~ temp.kT + consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Keeping only one explanatory variable: lme.ae.1a = lme(ae.logit ~ temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.1b = lme(ae.logit ~ consumer.type, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.1c = lme(ae.logit ~ log.weight, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Simply the mean: lme.ae.0 = lme(ae.logit ~ 1, data = ds.ae, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) ## Calculating the BIC and the delta-BIC for all models in the fixed effects selection part (Table S6): BIC(lme.ae.5,lme.ae.4a,lme.ae.4b,lme.ae.3a,lme.ae.3b,lme.ae.3c,lme.ae.2a,lme.ae.2b,lme.ae.2c,lme.ae.1a,lme.ae.1b,lme.ae.1c,lme.ae.0) BIC(lme.ae.5,lme.ae.4a,lme.ae.4b,lme.ae.3a,lme.ae.3b,lme.ae.3c,lme.ae.2a,lme.ae.2b,lme.ae.2c,lme.ae.1a,lme.ae.1b,lme.ae.1c,lme.ae.0)[2]-min(BIC(lme.ae.5,lme.ae.4a,lme.ae.4b,lme.ae.3a,lme.ae.3b,lme.ae.3c,lme.ae.2a,lme.ae.2b,lme.ae.2c,lme.ae.1a,lme.ae.1b,lme.ae.1c,lme.ae.0)[2]) # The model "lme.ae.2c" is the "best" according to BIC, with being more than 2 delta units apart from the next best one. We clearly decide on that. See methods in the main text for details. # The final model should be fitted by REML for a better estimation of p-values (see main text methods and references therein for details). # We moreover would like to directly visualize the values for the three constants (here intercepts at T0) and the common slope, therefor we add "-1" to the model formulation. lme.ae.fin = lme(ae.logit ~ consumer.type + temp.kT-1, data = ds.ae, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) ## The final result for assimilation efficiency (Table S9): summary(lme.ae.fin) ## Calculating the R-squared value according to "https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/" (last accessed on the 29.3.2017), see references therein and our methods for details. Note that we can not use the R-squared function to calculate the r-squared of the above used model with separated consumer type levels, we therefore have to use the standard way: lme.ae.fin.Rsquared = lme(ae.logit ~ consumer.type + temp.kT, data = ds.ae, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) rsquared(lme.ae.fin.Rsquared) ### Calculates the standard deviations for the consumer types: round(sd(resid(lme.ae.fin)[ds.ae$consumer.type == "carnivore"]),3) round(sd(resid(lme.ae.fin)[ds.ae$consumer.type == "detritivore"]),3) round(sd(resid(lme.ae.fin)[ds.ae$consumer.type == "herbivore"]),3) ### Calculates the confidence intervals for all parameters based on the final model: round(intervals(lme.ae.fin)$fixed,3) ### Calculates the confidence intervals for the intercepts based on the final model, but inverse logit transformed ### (epsion_0 = exp(epsilon_0_accent) / 1 + exp(epsilon_0_accent)): round(inv.logit(intervals(lme.ae.fin)$fixed[1:3]),3) ####################################################################################################################################### ################################################## Statistics for respiration rates ################################################### ####################################################################################################################################### ds.met = subset(ds.raw, metabolism.J.h > -999) # Select only data for assimilation efficiency (ds.met = dataset containing only assimilation efficiency). nrow(ds.met) # number of data points #### Selection for the best random stucture and heterogeneity - we use REML here, see Methods for details: lme.met.weights.REML = lme(log.met ~ (log.weight+temp.kT)*consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.REML = lme(log.met ~ (log.weight+temp.kT)*consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML" ) lme.met.weights.REML.b = lme(log.met ~ (log.weight+temp.kT)*consumer.type, data = ds.met, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.REML.b = lme(log.met ~ (log.weight+temp.kT)*consumer.type, data = ds.met, random = ~ 1 | reference.short, method = "REML" ) ## Calculating the BIC and the delta-BIC for the 4 REML models with varied random structure and data heterogeneity (Table S5): BIC(lme.met.weights.REML,lme.met.REML, lme.met.weights.REML.b,lme.met.REML.b) BIC(lme.met.weights.REML,lme.met.REML, lme.met.weights.REML.b,lme.met.REML.b)[2]-min(BIC(lme.met.weights.REML,lme.met.REML, lme.met.weights.REML.b,lme.met.REML.b)[2]) ## The model "lme.met.weights.REML.b" is the best, more than 2 delta-BIC units apart from the next best model - we choose this random structure and weights for further selection of the fixed effects. #### Selection for the best fixed structure, we use ML for model selection - see methods and references therein for details: # The full model: lme.met.5 = lme(log.met ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight + consumer.type:temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding one interaction between explanatory variables: lme.met.4a = lme(log.met ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.4b = lme(log.met ~ log.weight + temp.kT + consumer.type + consumer.type:temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, or one interaction and one explanatory variable: lme.met.3a = lme(log.met ~ log.weight + temp.kT + consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.3b = lme(log.met ~ log.weight + consumer.type + consumer.type:log.weight, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.3c = lme(log.met ~ temp.kT + consumer.type + consumer.type:temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, and one explanatory variable: lme.met.2a = lme(log.met ~ log.weight + temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.2b = lme(log.met ~ log.weight + consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.2c = lme(log.met ~ temp.kT + consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Keeping only one explanatory variable: lme.met.1a = lme(log.met ~ temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.1b = lme(log.met ~ consumer.type, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.met.1c = lme(log.met ~ log.weight, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Simply the mean: lme.met.0 = lme(log.met ~ 1, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) ## Calculating the delta-BIC for all models in the fixed effects selection part (Table S12): BIC(lme.met.5,lme.met.4a,lme.met.4b,lme.met.3a,lme.met.3b,lme.met.3c,lme.met.2a,lme.met.2b,lme.met.2c,lme.met.1a,lme.met.1b,lme.met.1c,lme.met.0) BIC(lme.met.5,lme.met.4a,lme.met.4b,lme.met.3a,lme.met.3b,lme.met.3c,lme.met.2a,lme.met.2b,lme.met.2c,lme.met.1a,lme.met.1b,lme.met.1c,lme.met.0)[2]-min(BIC(lme.met.5,lme.met.4a,lme.met.4b,lme.met.3a,lme.met.3b,lme.met.3c,lme.met.2a,lme.met.2b,lme.met.2c,lme.met.1a,lme.met.1b,lme.met.1c,lme.met.0)[2]) # The model "lme.met.2a" is the "best" according to BIC, with being more than 2 delta units apart from the next best one. We clearly decide on that. See methods in the main text for details. # The final model should be fitted by REML for a better estimation of p-values (see main text methods and references therein for details). lme.met.fin = lme(log.met ~ log.weight + temp.kT, data = ds.met, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) ## The final result for metabolic rates (Table S13): summary(lme.met.fin) ## The R-squared values according to "https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/" (last accessed on the 29.3.2017), see references therein and our methods for details: rsquared(lme.met.fin) ### Calculates the standard deviations for the consumer types: round(sd(resid(lme.met.fin)[ds.met$consumer.type == "carnivore"]),3) round(sd(resid(lme.met.fin)[ds.met$consumer.type == "detritivore"]),3) round(sd(resid(lme.met.fin)[ds.met$consumer.type == "herbivore"]),3) ### Calculates the confidence intervals for all parameters based on the final model: round(intervals(lme.met.fin)$fixed,3) ####################################################################################################################################### ################################################## Calculation of maintenance feeding rates ########################################### ####################################################################################################################################### #### First, we create an empty dataframe to save the results of the upcoming simulations in (CT=consumer type; C=carnivore; D = detritivore; H = herbivore): results.mf = data.frame(CTC = numeric(), CTD = numeric(), CTH = numeric(), CTC.log.weight = numeric(), CTD.log.weight = numeric(), CTH.log.weight = numeric(), CTC.temp.kT = numeric(), CTD.temp.kT = numeric(), CTH.temp.kT = numeric(), sim = integer() ) #### Second, we create an empty dataframe to save the simulated data in: ds.mf = data.frame(log.weight = numeric(), temp.kT = numeric(), consumer.type = factor(), log.mf = numeric(), sim = integer() ) set.seed(667) # we set a seed to create reproducible results for(i in 1:1000){ ########################################################################################################################## ######################################################## Step (1) ######################################################## ###################################### Select number of data points to be simulated ###################################### ########################################################################################################################## #### In the following we sample three times the number of datapoints we want to simulate below. #### The first number is for consumer (no.C), the second for herbivores (no.H), and the third for detritivores (no.D). #### We use the number of datapoints we had in the original data. The data of assimilation was always less than the data from the respiration data. #### Therefore, the lower end from where to sample is always from the assimilation data, the upper from the respiration data. no.C = sample(nrow(ds.ae[ds.ae$consumer.type=="carnivore",]):nrow(ds.met[ds.met$consumer.type=="carnivore",]),1) no.D = sample(nrow(ds.ae[ds.ae$consumer.type=="detritivore",]):nrow(ds.met[ds.met$consumer.type=="detritivore",]),1) no.H = sample(nrow(ds.ae[ds.ae$consumer.type=="herbivore",]):nrow(ds.met[ds.met$consumer.type=="herbivore",]),1) ########################################################################################################################## ######################################################## Step (2) ######################################################## ########################################### Sample body masses and temperatures ########################################## ########################################################################################################################## #### Here we sample Arrhenius temperatures out of the original datasets; the number of datapoints per consumer type equals the number sampled above: temp.kT.C = sample(ds.raw$temp.kT[ds.raw$consumer.type=="carnivore"], no.C) temp.kT.D = sample(ds.raw$temp.kT[ds.raw$consumer.type=="detritivore"], no.D) temp.kT.H = sample(ds.raw$temp.kT[ds.raw$consumer.type=="herbivore"], no.H) #### Here we sample log.weights out of the original datasets; the number of datapoints per consumer type equals the number sampled above: mass.C = sample(ds.raw$log.weight[ds.raw$consumer.type=="carnivore"], no.C) mass.D = sample(ds.raw$log.weight[ds.raw$consumer.type=="detritivore"], no.D) mass.H = sample(ds.raw$log.weight[ds.raw$consumer.type=="herbivore"], no.H) #### Create a integrated dataset to analyze (mf = maintenance feeding): ds.mf.i = data.frame(log.weight = c(mass.C,mass.D,mass.H), temp.kT = c(temp.kT.C,temp.kT.D,temp.kT.H), consumer.type = c(rep("carnivore",no.C),rep("detritivore",no.D),rep("herbivore",no.H)) ) ########################################################################################################################## ######################################################## Step (3) ######################################################## ############################### Sample intercepts and slopes from the statistical results ################################ ########################################################################################################################## ## We select an intercept for each consumer type, according to the estimate of the "best" assimilation efficiency fitting, varied by the standard error of the statistical results: sim.i.C.ae = rnorm(1, mean = summary(lme.ae.fin)$tTable[1], sd = summary(lme.ae.fin)$tTable[5]) sim.i.D.ae = rnorm(1, mean = summary(lme.ae.fin)$tTable[2], sd = summary(lme.ae.fin)$tTable[6]) sim.i.H.ae = rnorm(1, mean = summary(lme.ae.fin)$tTable[3], sd = summary(lme.ae.fin)$tTable[7]) ## We select a shared slope according to the assimilation efficiency fitting, varied by the standard error of the statistical result: sim.t.ae = rnorm(1, mean = summary(lme.ae.fin)$tTable[4], sd = summary(lme.ae.fin)$tTable[8]) ## We now repeat the same procedure for respiration rates; we calculate a shared intercept, allometric slope and activation energy: sim.i.met = rnorm(1, mean = summary(lme.met.fin)$tTable[1], sd = summary(lme.met.fin)$tTable[4]) sim.m.met = rnorm(1, mean = summary(lme.met.fin)$tTable[2], sd = summary(lme.met.fin)$tTable[5]) sim.t.met = rnorm(1, mean = summary(lme.met.fin)$tTable[3], sd = summary(lme.met.fin)$tTable[6]) ########################################################################################################################## ######################################################## Step (4) ######################################################## ################################## Calculate assimilation efficiency and metabolic rates ################################# ########################################################################################################################## ## We calculate the assimilation efficiencies for the sampled temperatues and weights (no effect here) for all three consumer types. ## We additionally add variation to the data using the standard deviation of the conditional residuals for each consumer type: sim.C.ae = sim.i.C.ae + sim.t.ae*temp.kT.C + rnorm(no.C, mean = 0, sd = sd(as.vector(resid(lme.ae.fin))[ds.ae$consumer.type == "carnivore"])) sim.D.ae = sim.i.D.ae + sim.t.ae*temp.kT.D + rnorm(no.D, mean = 0, sd = sd(as.vector(resid(lme.ae.fin))[ds.ae$consumer.type == "detritivore"])) sim.H.ae = sim.i.H.ae + sim.t.ae*temp.kT.H + rnorm(no.H, mean = 0, sd = sd(as.vector(resid(lme.ae.fin))[ds.ae$consumer.type == "herbivore"])) ## We calculate the respiration rates for the sampled temperatues and weights for all three consumer types (no effect here). ## We additionally add variation to the data using the standard deviation of the conditional residuals for each consumer type: sim.C.met = sim.i.met + sim.t.met*temp.kT.C + sim.m.met*mass.C + rnorm(no.C, mean = 0, sd= sd(as.vector(resid(lme.met.fin))[ds.met$consumer.type == "carnivore"])) sim.D.met = sim.i.met + sim.t.met*temp.kT.D + sim.m.met*mass.D + rnorm(no.D, mean = 0, sd= sd(as.vector(resid(lme.met.fin))[ds.met$consumer.type == "detritivore"])) sim.H.met = sim.i.met + sim.t.met*temp.kT.H + sim.m.met*mass.H + rnorm(no.H, mean = 0, sd= sd(as.vector(resid(lme.met.fin))[ds.met$consumer.type == "herbivore"])) ########################################################################################################################## ######################################################## Step (5) ######################################################## ############################################## Calculate maintenance feeding ############################################# ########################################################################################################################## ## Here we calculate the maintenance feeding (mf) and add them to the dataset "ds.mf.i": ds.mf.i$log.mf = c(log(exp(sim.C.met)/gtools::inv.logit(sim.C.ae)), log(exp(sim.D.met)/gtools::inv.logit(sim.D.ae)), log(exp(sim.H.met)/gtools::inv.logit(sim.H.ae))) ########################################################################################################################## ######################################################## Step (6) ######################################################## ####################################### Statistical analyses of maintenance feeding ###################################### ########################################################################################################################## #### Estimating the intercepts, allometric slopes and activation energies using a simple lm, we don't have random effects here: lm.mf = lm(log.mf ~ consumer.type/(log.weight + temp.kT) -1, data = ds.mf.i) ## for easier intercept and slope extracting! ## Saving the results in a dataframe that is only created in the for-loop: results.mf.i = data.frame(CTC = coef(lm.mf)[[1]], CTD = coef(lm.mf)[[2]], CTH = coef(lm.mf)[[3]], CTC.log.weight = coef(lm.mf)[[4]], CTD.log.weight = coef(lm.mf)[[5]], CTH.log.weight = coef(lm.mf)[[6]], CTC.temp.kT = coef(lm.mf)[[7]], CTD.temp.kT = coef(lm.mf)[[8]], CTH.temp.kT = coef(lm.mf)[[9]], sim = i ) ## Next, we bind the results-dataframe to the overall dataframe: results.mf = rbind(results.mf,results.mf.i) ## We add the information which simulation the data belongs to to the data-dataframe: ds.mf.i$sim = rep(i,nrow(ds.mf.i)) ## Subsequently, we bind the internal for-loop data to the overall dataset: ds.mf = rbind(ds.mf,ds.mf.i) print(i) ## Prints out where we are in the loop. } ## we re-organize the dataframe results.mf for upcoming statistics: results.mf2 = data.frame(i = c(results.mf$CTC , results.mf$CTD , results.mf$CTH), a = c(results.mf$CTC.log.weight, results.mf$CTD.log.weight, results.mf$CTH.log.weight), E = c(results.mf$CTC.temp.kT , results.mf$CTD.temp.kT , results.mf$CTH.temp.kT), CT = c(rep("C",1000),rep("D",1000),rep("H",1000)) ) str(results.mf2) ### Here we test if the intercepts are shared across consumer types or if all share the same intercept: mod.i1 = lm(i ~ CT, data = results.mf2, na.action = "na.fail") mod.i2 = lm(i ~ 1, data = results.mf2, na.action = "na.fail") BIC(mod.i1,mod.i2)[2]-min(BIC(mod.i1,mod.i2)[2]) ## clear signal!! mod.i.fin = lm(i ~ CT-1, data = results.mf2, na.action = "na.fail") summary(mod.i.fin) ### Here we test if the activation energy is shared across consumer types or if all share the same slope mod.E1 = lm(E ~ CT, data = results.mf2, na.action = "na.fail") mod.E2 = lm(E ~ 1, data = results.mf2, na.action = "na.fail") BIC(mod.E1,mod.E2)[2]-min(BIC(mod.E1,mod.E2)[2]) mod.E.fin = lm(E ~ CT-1, data = results.mf2, na.action = "na.fail") summary(mod.E.fin) ### Here we test if the allometric slope is shared across consumer types or if all share the same slope: mod.a1 = lm(a ~ CT, data = results.mf2, na.action = "na.fail") mod.a2 = lm(a ~ 1, data = results.mf2, na.action = "na.fail") BIC(mod.a1,mod.a2)[2]-min(BIC(mod.a1,mod.a2)[2]) ## the dBIC reveals that we only need a shared allometric slope summary(mod.a2) ####################################################################################################################################### ################################################## Fixed Effects Statistics for assimilation efficiency ############################### ################################################## using the alternative random structure ############################################# ####################################################################################################################################### #### Selection for the best fixed structure, we use ML for model selection - see methods and references therein for details: # The full model: lme.ae.5.sensitivity = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding one interaction between explanatory variables: lme.ae.4a.sensitivity = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.4b.sensitivity = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, or one interaction and one explanatory variable: lme.ae.3a.sensitivity = lme(ae.logit ~ log.weight + temp.kT + consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.3b.sensitivity = lme(ae.logit ~ log.weight + consumer.type + consumer.type:log.weight, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.3c.sensitivity = lme(ae.logit ~ temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, and one explanatory variable: lme.ae.2a.sensitivity = lme(ae.logit ~ log.weight + temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.2b.sensitivity = lme(ae.logit ~ log.weight + consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.2c.sensitivity = lme(ae.logit ~ temp.kT + consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Keeping only one explanatory variable: lme.ae.1a.sensitivity = lme(ae.logit ~ temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.1b.sensitivity = lme(ae.logit ~ consumer.type, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.1c.sensitivity = lme(ae.logit ~ log.weight, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Simply the mean: lme.ae.0.sensitivity = lme(ae.logit ~ 1, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) ## Calculating the BIC and the delta-BIC for all models in the fixed effects selection part (Table S7): BIC(lme.ae.5.sensitivity,lme.ae.4a.sensitivity,lme.ae.4b.sensitivity,lme.ae.3a.sensitivity,lme.ae.3b.sensitivity,lme.ae.3c.sensitivity,lme.ae.2a.sensitivity,lme.ae.2b.sensitivity,lme.ae.2c.sensitivity,lme.ae.1a.sensitivity,lme.ae.1b.sensitivity,lme.ae.1c.sensitivity,lme.ae.0.sensitivity) BIC(lme.ae.5.sensitivity,lme.ae.4a.sensitivity,lme.ae.4b.sensitivity,lme.ae.3a.sensitivity,lme.ae.3b.sensitivity,lme.ae.3c.sensitivity,lme.ae.2a.sensitivity,lme.ae.2b.sensitivity,lme.ae.2c.sensitivity,lme.ae.1a.sensitivity,lme.ae.1b.sensitivity,lme.ae.1c.sensitivity,lme.ae.0.sensitivity)[2]-min(BIC(lme.ae.5.sensitivity,lme.ae.4a.sensitivity,lme.ae.4b.sensitivity,lme.ae.3a.sensitivity,lme.ae.3b.sensitivity,lme.ae.3c.sensitivity,lme.ae.2a.sensitivity,lme.ae.2b.sensitivity,lme.ae.2c.sensitivity,lme.ae.1a.sensitivity,lme.ae.1b.sensitivity,lme.ae.1c.sensitivity,lme.ae.0.sensitivity)[2]) # The model "lme.ae.2c" is the "best" according to BIC, with being more than 2 delta units apart from the next best one. We clearly decide on that. See methods in the main text for details. # The final model should be fitted by REML for a better estimation of p-values (see main text methods and references therein for details). # We moreover would like to directly visualize the values for the three constants (here intercepts at T0) and the common slope, therefor we add "-1" to the model formulation. lme.ae.fin.sensitivity = lme(ae.logit ~ consumer.type + temp.kT-1, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) ## The final result for assimilation efficiency (Table S10): summary(lme.ae.fin.sensitivity) ## Calculating the R-squared value according to "https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/" (last accessed on the 29.3.2017), see references therein and our methods for details. Note that we can not use the R-squared function to calculate the r-squared of the above used model with separated consumer type levels, we therefore have to use the standard way: lme.ae.fin.sensitivity.Rsquared = lme(ae.logit ~ consumer.type + temp.kT, data = ds.ae, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) rsquared(lme.ae.fin.sensitivity.Rsquared) ### Calculates the standard deviations for the consumer types: round(sd(resid(lme.ae.fin.sensitivity)[ds.ae$consumer.type == "carnivore"]),3) round(sd(resid(lme.ae.fin.sensitivity)[ds.ae$consumer.type == "detritivore"]),3) round(sd(resid(lme.ae.fin.sensitivity)[ds.ae$consumer.type == "herbivore"]),3) ### Calculates the confidence intervals for all parameters based on the final model: round(intervals(lme.ae.fin.sensitivity)$fixed,3) ### Calculates the confidence intervals for the intercepts based on the final model, but inverse logit transformed ### (epsion_0 = exp(epsilon_0_accent) / 1 + exp(epsilon_0_accent)): round(inv.logit(intervals(lme.ae.fin.sensitivity)$fixed[1:3]),3) ####################################################################################################################################### ################################################## Fixed Effects Statistics for assimilation efficiency ############################### ################################################## using only studies where fresh weights were reported ############################### ####################################################################################################################################### ds.ae.fw = subset(ds.ae, fresh.weight.calculation == "original") # Select only data from studies reported fresh weight case for assimilation efficiency (ds.ae.fw = dataset containing only assimilation efficiency from studies reporting fresh weights). nrow(ds.ae.fw) # number of data points #### Selection for the best random stucture and heterogeneity - we use REML here, see Methods for details: lme.ae.fw.weights.REML = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae.fw, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.REML = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae.fw, random = ~ 1 | taxonomic.group.consumer/reference.short, method = "REML" ) lme.ae.fw.weights.REML.b = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.REML.b = lme(ae.logit ~ (log.weight+temp.kT)*consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "REML" ) ## Calculating the delta-BIC for the 4 REML models with varied random structure and data heterogeneity (Table S4): BIC(lme.ae.fw.weights.REML,lme.ae.fw.REML, lme.ae.fw.weights.REML.b,lme.ae.fw.REML.b) BIC(lme.ae.fw.weights.REML,lme.ae.fw.REML, lme.ae.fw.weights.REML.b,lme.ae.fw.REML.b)[2]-min(BIC(lme.ae.fw.weights.REML,lme.ae.fw.REML, lme.ae.fw.weights.REML.b,lme.ae.fw.REML.b)[2]) ## The model "lme.ae.fw.weights.REML.b" is the best, more than 2 delta-BIC units apart from the next best model - we choose this random structure with weights for further selection of the fixed effects. #### Selection for the best fixed structure, we use ML for model selection - see methods and references therein for details: # The full model: lme.ae.fw.5 = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight + consumer.type:temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding one interaction between explanatory variables: lme.ae.fw.4a = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:log.weight, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.4b = lme(ae.logit ~ log.weight + temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, or one interaction and one explanatory variable: lme.ae.fw.3a = lme(ae.logit ~ log.weight + temp.kT + consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.3b = lme(ae.logit ~ log.weight + consumer.type + consumer.type:log.weight, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.3c = lme(ae.logit ~ temp.kT + consumer.type + consumer.type:temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Excluding both interactions between explanatory variables, and one explanatory variable: lme.ae.fw.2a = lme(ae.logit ~ log.weight + temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.2b = lme(ae.logit ~ log.weight + consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.2c = lme(ae.logit ~ temp.kT + consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Keeping only one explanatory variable: lme.ae.fw.1a = lme(ae.logit ~ temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.1b = lme(ae.logit ~ consumer.type, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) lme.ae.fw.1c = lme(ae.logit ~ log.weight, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) # Simply the mean: lme.ae.fw.0 = lme(ae.logit ~ 1, data = ds.ae.fw, random = ~ 1 | reference.short, method = "ML", weights = varIdent(form = ~1 | consumer.type) ) ## Calculating the BIC and the delta-BIC for all models in the fixed effects selection part (Table S8): BIC(lme.ae.fw.5,lme.ae.fw.4a,lme.ae.fw.4b,lme.ae.fw.3a,lme.ae.fw.3b,lme.ae.fw.3c,lme.ae.fw.2a,lme.ae.fw.2b,lme.ae.fw.2c,lme.ae.fw.1a,lme.ae.fw.1b,lme.ae.fw.1c,lme.ae.fw.0) BIC(lme.ae.fw.5,lme.ae.fw.4a,lme.ae.fw.4b,lme.ae.fw.3a,lme.ae.fw.3b,lme.ae.fw.3c,lme.ae.fw.2a,lme.ae.fw.2b,lme.ae.fw.2c,lme.ae.fw.1a,lme.ae.fw.1b,lme.ae.fw.1c,lme.ae.fw.0)[2]-min(BIC(lme.ae.fw.5,lme.ae.fw.4a,lme.ae.fw.4b,lme.ae.fw.3a,lme.ae.fw.3b,lme.ae.fw.3c,lme.ae.fw.2a,lme.ae.fw.2b,lme.ae.fw.2c,lme.ae.fw.1a,lme.ae.fw.1b,lme.ae.fw.1c,lme.ae.fw.0)[2]) # The model "lme.ae.fw.2c" is the "best" according to BIC, with being more than 2 delta units apart from the next best one. We clearly decide on that. See methods in the main text for details. # The final model should be fitted by REML for a better estimation of p-values (see main text methods and references therein for details). # We moreover would like to directly visualize the values for the three constants (here intercepts at T0) and the common slope, therefor we add "-1" to the model formulation. lme.ae.fw.fin = lme(ae.logit ~ consumer.type + temp.kT-1, data = ds.ae.fw, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) ## The final result for assimilation efficiency (Table S11): summary(lme.ae.fw.fin) ## Calculating the R-squared value according to "https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/" (last accessed on the 29.3.2017), see references therein and our methods for details. Note that we can not use the R-squared function to calculate the r-squared of the above used model with separated consumer type levels, we therefore have to use the standard way: lme.ae.fw.fin.Rsquared = lme(ae.logit ~ consumer.type + temp.kT, data = ds.ae.fw, random = ~ 1 | reference.short, method = "REML", weights = varIdent(form = ~1 | consumer.type) ) rsquared(lme.ae.fw.fin.Rsquared) ### Calculates the standard deviations for the consumer types: round(sd(resid(lme.ae.fw.fin)[ds.ae.fw$consumer.type == "carnivore"]),3) round(sd(resid(lme.ae.fw.fin)[ds.ae.fw$consumer.type == "detritivore"]),3) round(sd(resid(lme.ae.fw.fin)[ds.ae.fw$consumer.type == "herbivore"]),3) ### Calculates the confidence intervals for all parameters based on the final model: round(intervals(lme.ae.fw.fin)$fixed,3) ### Calculates the confidence intervals for the intercepts based on the final model, but inverse logit transformed ### (epsion_0 = exp(epsilon_0_accent) / 1 + exp(epsilon_0_accent)): round(inv.logit(intervals(lme.ae.fw.fin)$fixed[1:3]),3)