rm(list=ls(all=T)) # clears workspace set.seed(7) # Set library(metafor) library(brms) library(tidybayes) library(ggplot2) library(tidyr) ## Meta-analysis: effect of extreme events ## JD Gonzalez-Trujillo, APR2022 ## # 1. Database loading # 2. Meta-analysis using Hedges q # 3. Meta-analysis using Log-Ratio (only effect size estimation, repeat code in 2 for running the models) # 1. Database loading ----------------------------------------------------- #Loading databases # Estructural responses str_resp <- read.table("structural_resp.txt",h=T,dec=",",sep="\t", fileEncoding="latin1", stringsAsFactors=T) str_resp$year <- as.numeric(gsub("^.*([0-9]{4}).*", "\\1", str_resp$authors)) str(str_resp) # 198 effect sizes # Functional responses fun_resp <- read.table("functional_resp.txt",h=T,dec=",",sep="\t", fileEncoding="latin1", stringsAsFactors=T) fun_resp$year <- as.numeric(gsub("^.*([0-9]{4}).*", "\\1", fun_resp$authors)) str(fun_resp) # 134 effect sizes # Data for mapping the studies map_df <- read.csv2("refMap.csv",h=T,sep=",",dec=".", check.names = FALSE) head(map_df) # 2. Meta-analysis using Hedges g -------------------------------------------- # function to compute the variance corrected by sample size, # based on treatment means, standard deviations, and sample sizes # vi_n_es=function(mean_C,mean_E,n_C,n_E,sd_C,sd_E){ s=sqrt(((n_C-1)*(sd_C^2)+(n_E-1)*(sd_E^2))/(n_C+n_E-2)) m=(n_C+n_E-2) J=1-(3/(4*m-1)) n_tilde=n_E*n_C/(n_E+n_C) d=((mean_E-mean_C)/s)*J var_d_n=((1-3/(4*(n_E+n_C-2)-1))^2)*(n_E+n_C-2)/(n_tilde*(n_E+n_C-4)) return(c(var_d_n)) } #Estimating Hedges g - Ecosystem structure es_str <- escalc(measure = "SMD",vtype="UB", m1i = treat_avg,sd1i = treat_sd,n1i = treat_n, m2i = ctrl_avg, sd2i = ctrl_sd, n2i = ctrl_n,data = str_resp, slab=str_resp$authors) # variance corrected by sample size es_str$vi_n <- vi_n_es(mean_C = str_resp$ctrl_avg,mean_E = str_resp$treat_avg, n_C = str_resp$ctrl_n,n_E = str_resp$treat_n, sd_C = str_resp$ctrl_sd,sd_E = str_resp$treat_sd) #Estimating Hedges g - ecosystem functioning es_fun <- escalc(measure = "SMD",vtype="UB", m1i = treat_avg,sd1i = treat_sd,n1i = treat_n, m2i = ctrl_avg, sd2i = ctrl_sd, n2i = ctrl_n,data = fun_resp, slab=fun_resp$authors) # variance corrected by sample size es_fun$vi_n <- vi_n_es(mean_C = fun_resp$ctrl_avg,mean_E = fun_resp$treat_avg, n_C = fun_resp$ctrl_n,n_E = fun_resp$treat_n, sd_C = fun_resp$ctrl_sd,sd_E = fun_resp$treat_sd) # 2a. Species richness str_S <- subset(es_str,es_str$variable == "Richness") table(str_S$biological.group) # retaining biological groups with 3 or more studies #str_S <- subset(str_S,str_S$biological.group == "invert") # keeping invertebrates str_S <- droplevels(str_S) str(str_S) # Number of efect sizes table(str_S$Climate.anomaly) # 41 effect sizes # Create a funnel plot and perform Egger's test # funnel(x=str_S$yi, sei=str_S$vi, pch=15,main="A. Species richness") regtest(x= str_S$yi, sei= str_S$vi_n) # Example: Find the best random effect structure for the Hedges' g data (Repeat for every model below) # Study_ID b0 <- brm(data = str_S, family = gaussian, yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1 | study_ID), iter = 10000, warmup = 1000, cores = 3, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99), prior = c(prior(normal(0, 10), "b"), prior(cauchy(0, 1), "sd"))) # Study duration b1 <- brm(data = str_S, family = gaussian, yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1 | Duration), iter = 10000, warmup = 1000, cores = 3, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99), prior = c(prior(normal(0, 10), "b"), prior(cauchy(0, 1), "sd"))) # Study type b2 <- brm(data = str_S, family = gaussian, yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1| Study.type), iter = 10000, warmup = 1000, cores = 3, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99), prior = c(prior(normal(0, 10), "b"), prior(cauchy(0, 1), "sd"))) # Study ID + type b3 <- brm(data = str_S, family = gaussian, yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1| Study.type)+ (1 | study_ID), iter = 10000, warmup = 1000, cores =3, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99), prior = c(prior(normal(0, 10), "b"))) # Do leave-one-out cross-validation for model selection b0_loo<-loo(b0, reloo=T)# The best b1_loo<-loo(b1, reloo=T) b2_loo<-loo(b2, reloo=T) b3_loo<-loo(b3, reloo=T) # Compare the LOOIC print(loo_compare(b0_loo, b1_loo,b2_loo,b3_loo), simplify=FALSE) # model b1 with random effects of Study_ID is best supported # mixed models ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_S <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = str_S, family = gaussian(), iter = 10000, warmup = 1000, cores = 4, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_S) bm_S_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1 + (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = str_S, family = gaussian(), iter = 10000, warmup = 1000, cores = 4, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_S_t) # 2b. Biomass str_bioM <- subset(es_str,es_str$variable == "Biomass") table(str_bioM$biological.group) # retaining biological groups with 3 or more studies str_bioM <- subset(str_bioM,str_bioM$biological.group != "bacteria") str_bioM <- subset(str_bioM,str_bioM$biological.group != "Protozoa") str_bioM <- subset(str_bioM,str_bioM$biological.group != "fish") #str_bioM <- subset(str_bioM,str_bioM$biological.group == "invert") str_bioM <- droplevels(str_bioM) table(str_bioM$Climate.anomaly) # Create a funnel plot and perform Egger's test # str_bioM <- str_bioM[str_bioM$yi > c(-9),] funnel(x=str_bioM$yi, sei=str_bioM$vi,main="B. Biomass") regtest(x= str_bioM$yi, sei= str_bioM$vi_n) table(str_bioM$Climate.anomaly) # mixed models ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_bioM <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = str_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM) bm_bioM_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = str_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_t) # mixed models - Biofilm biomass biof_bioM <- subset(str_bioM,str_bioM$biological.group == "biofilm") ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_bioM_biof <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = biof_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_biof) bm_bioM_biof_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = biof_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_biof_t) # mixed models - Fungi biomass fungi_bioM <- subset(str_bioM,str_bioM$biological.group == "fungi") ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_bioM_fungi <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = fungi_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_fungi) bm_bioM_fungi_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = fungi_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_fungi_t) # mixed models - Invertebrate biomass inv_bioM <- subset(str_bioM,str_bioM$biological.group == "invert") ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_bioM_inv <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = inv_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_inv) bm_bioM_inv_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1+ (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = inv_bioM, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_bioM_inv_t) # 2c. Functional - decomposition fun_db <- es_fun table(fun_db$function.) # retaining biological groups with 3 or more studies fun_decom <- subset(fun_db,fun_db$function. == "decomposition") # removing growth table(fun_decom$Climate.anomaly) # Create a funnel plot and perform Egger's test # fun_decom <- fun_decom[fun_decom$yi < c(10),] funnel(x=fun_decom$yi, sei=fun_decom$vi, pch=17,main="D. Decomposition") regtest(x= fun_decom$yi, sei= fun_decom$vi_n) # model b1 with random effects of Study_ID is best supported # mixed models ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_decom <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = fun_decom, family = gaussian(), iter = 10000, warmup = 1000, cores = 4, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_decom) bm_decom_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1 + (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = fun_decom, family = gaussian(), iter = 10000, warmup = 1000, cores = 4, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_decom_t) # 2d. Functional - respiration fun_metab <- subset(fun_db,fun_db$function. == "Respiration") # removing growth table(fun_metab$Climate.anomaly) # Create a funnel plot and perform Egger's test # fun_metab <- fun_metab[fun_metab$yi > c(-9),] funnel(x=fun_metab$yi, sei=fun_metab$vi, pch=15, main = "E. Respiration") regtest(x= fun_metab$yi, sei= fun_metab$vi_n) # mixed models ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_metab <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1+ (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = fun_metab, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_metab) bm_metab_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1+ (1 | study_ID), prior = c(prior(normal(0, 10), "b"), prior(cauchy( 0, 1), "sd")), data = fun_metab, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_metab_t) # 2d. Functional - Primary productivity fun_metab <- subset(fun_db,fun_db$function. == "Primary Productivity") # removing growth # Create a funnel plot and perform Egger's test # funnel(x=fun_metab$yi, sei=fun_metab$vi, pch=15, main="F. Primary productivity") regtest(x= fun_metab$yi, sei= fun_metab$vi_n) # mixed models ## weighted using the sample size rather than within-study variance ## as recommended by Hamman et al. 2018. Ecosphere 9(9):e02419. 10.1002/ecs2.2419 bm_pp <- brm(yi | se(vi_n, sigma = TRUE) ~ Climate.anomaly - 1 + (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = fun_metab, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_pp) bm_pp_t <- brm(yi | se(vi_n, sigma = TRUE) ~ Anomaly.type - 1 + (1| study_ID), prior = c(prior(normal(0, 10), "b")), data = fun_metab, family = gaussian(), iter = 10000, warmup = 1000, cores = 6, chains = 4, seed = 14, control = list(max_treedepth = 12, adapt_delta = 0.99)) summary(bm_pp_t) # 3. Log response ratio computation -------------------------------------------- #Estimating Log-response ratio - Ecosystem structure es_str_lrr <- escalc(measure = "ROM", m1i = treat_avg,sd1i = treat_sd,n1i = treat_n, m2i = ctrl_avg, sd2i = ctrl_sd, n2i = ctrl_n,data = str_resp, slab=str_resp$authors,var.names=c("LRR","LRR_var"),digits=4) #Estimating Log-response ratio - ecosystem functioning es_fun_lrr <- escalc(measure = "ROM", m1i = treat_avg,sd1i = treat_sd,n1i = treat_n, m2i = ctrl_avg, sd2i = ctrl_sd, n2i = ctrl_n,data = fun_resp, slab=fun_resp$authors,var.names=c("LRR","LRR_var"),digits=4) # 2a. Species richness str_S_lrr <- subset(es_str_lrr,es_str$variable == "Richness") table(str_S$biological.group) # retaining biological groups with 3 or more studies #str_S <- subset(str_S,str_S$biological.group == "invert") # keeping invertebrates str_S <- droplevels(str_S) str(str_S)