Rerun of Clements et al data corrected and complete 0.1 and 0.001 remove low CO2 under 800
Anon
18/07/2022
—————————– Supplementary information to —————————-
—————— META-ANALYSIS REVEALS AN EXTREME “DECLINE EFFECT” —————— —————— IN OCEAN ACIDIFICATION IMPACTS ON FISH BEHAVIOUR ——————-
——– Jeff C. Clements, Josefin Sundin, Timothy D. Clark, Fredrik Jutfelt ———-
Get packages
library(pacman)## Warning: package 'pacman' was built under R version 4.1.2pacman::p_load(metafor, MCMCglmm, tidyverse, rotl, magrittr, kableExtra, rmarkdown,gridExtra, psych, bindrcpp, pander)
library(BiocManager)## Warning: package 'BiocManager' was built under R version 4.1.2## Bioconductor version '3.14' is out-of-date; the current release version '3.15'
## is available with R version '4.2'; see https://bioconductor.org/installlibrary(ggplot2)
library(viridis)## Loading required package: viridisLitelibrary(patchwork)## Warning: package 'patchwork' was built under R version 4.1.2META-ANALYSIS - YEAR ONLINE - FULL DATASET
##attach dataset
decline<-read.csv(file.choose()) ##use dataset "S5 Data"
attach(decline)##set factors
decline$year.online<-as.factor(decline$year.online)
decline$year.print<-as.factor(decline$year.print)
decline$obs<-as.factor(decline$obs)
decline$study<-as.factor(decline$study)##view summary
summary(decline)## obs study authors year.online year.print
## 1 : 1 a87 : 40 Length:644 2018 :148 2018 :168
## 2 : 1 a90 : 36 Class :character 2015 : 99 2015 : 85
## 3 : 1 a31 : 28 Mode :character 2017 : 72 2016 : 71
## 4 : 1 a22 : 24 2014 : 65 2013 : 66
## 5 : 1 a73 : 22 2013 : 54 2014 : 58
## 6 : 1 a42 : 21 2019 : 52 2017 : 57
## (Other):638 (Other):473 (Other):154 (Other):139
## if.at.pub X2017.if if.group avg.n
## Length:644 Length:644 Length:644 Min. : 4.25
## Class :character Class :character Class :character 1st Qu.: 11.95
## Mode :character Mode :character Mode :character Median : 18.10
## Mean : 26.69
## 3rd Qu.: 30.00
## Max. :156.17
## NA's :553
## species climate cue cue.type
## Length:644 Length:644 Length:644 Length:644
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## life.stage ctrl.n ctrl.mean ctrl.sd
## Length:644 Min. : 3.00 Min. : -61.58 Min. : 0.000
## Class :character 1st Qu.: 10.00 1st Qu.: 1.39 1st Qu.: 1.197
## Mode :character Median : 18.00 Median : 10.32 Median : 5.867
## Mean : 27.73 Mean : 548.23 Mean : 127.582
## 3rd Qu.: 36.00 3rd Qu.: 46.10 3rd Qu.: 21.804
## Max. :245.00 Max. :154936.88 Max. :25490.446
##
## oa.n oa.mean oa.sd
## Min. : 2.00 Min. : -25.51 Min. : 0.00
## 1st Qu.: 10.00 1st Qu.: 1.44 1st Qu.: 1.08
## Median : 18.00 Median : 13.19 Median : 6.47
## Mean : 27.64 Mean : 546.48 Mean : 138.96
## 3rd Qu.: 36.00 3rd Qu.: 45.52 3rd Qu.: 20.33
## Max. :249.00 Max. :157061.25 Max. :36812.37
## ##subset by year
y2009 <- filter(decline, year.online == "2009")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2009$obs <- 1:nrow(y2009)
y2010 <- filter(decline, year.online == "2010")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2010$obs <- 1:nrow(y2010)
y2011 <- filter(decline, year.online == "2011")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2011$obs <- 1:nrow(y2011)
y2012 <- filter(decline, year.online == "2012")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2012$obs <- 1:nrow(y2012)
y2013 <- filter(decline, year.online == "2013")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2013$obs <- 1:nrow(y2013)
y2014 <- filter(decline, year.online == "2014")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2014$obs <- 1:nrow(y2014)
y2015 <- filter(decline, year.online == "2015")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2015$obs <- 1:nrow(y2015)
y2016 <- filter(decline, year.online == "2016")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2016$obs <- 1:nrow(y2016)
y2017 <- filter(decline, year.online == "2017")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2017$obs <- 1:nrow(y2017)
y2018 <- filter(decline, year.online == "2018")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2018$obs <- 1:nrow(y2018)
y2019 <- filter(decline, year.online == "2019")[,-match(c("avg.n","cue.type","if.at.pub","X2017.if","if.group"), colnames (decline))]
y2019$obs <- 1:nrow(y2019)##compute effect sizes for each year
lnRR2009 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2009,append=TRUE)
lnRR2010 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2010,append=TRUE)
lnRR2011 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2011,append=TRUE)
lnRR2012 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2012,append=TRUE)
lnRR2013 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2013,append=TRUE)
lnRR2014 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2014,append=TRUE)
lnRR2015 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2015,append=TRUE)
lnRR2016 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2016,append=TRUE)
lnRR2017 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2017,append=TRUE)
lnRR2018 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2018,append=TRUE)
lnRR2019 <- escalc(measure= "ROM", m1i=oa.mean,sd1i=oa.sd,n1i=oa.n,m2i=ctrl.mean,sd2i=ctrl.sd,n2i=ctrl.n,data=y2019,append=TRUE)
#Note log(m1i/m21i) produced NAs for 2011, 2012, 2013, 2015, 2017##remove NAs
lnRR2009clean<-na.omit(lnRR2009)
lnRR2010clean<-na.omit(lnRR2010)
lnRR2011clean<-na.omit(lnRR2011)
lnRR2012clean<-na.omit(lnRR2012)
lnRR2013clean<-na.omit(lnRR2013)
lnRR2014clean<-na.omit(lnRR2014)
lnRR2015clean<-na.omit(lnRR2015)
lnRR2016clean<-na.omit(lnRR2016)
lnRR2017clean<-na.omit(lnRR2017)
lnRR2018clean<-na.omit(lnRR2018)
lnRR2019clean<-na.omit(lnRR2019)##view mean-variance relationship
pp1<-ggplot(decline,aes(x=log(ctrl.mean),y=log(ctrl.sd),col=year.online))+ geom_point(size=2,na.rm = TRUE)
pp2<-ggplot(decline,aes(x=log(oa.mean),y=log(oa.sd),col=year.online))+ geom_point(size=2,na.rm = TRUE)
grid.arrange(pp1,pp2, nrow =1)## Warning in log(ctrl.mean): NaNs produced
## Warning in log(ctrl.mean): NaNs produced## Warning in log(oa.mean): NaNs produced
## Warning in log(oa.mean): NaNs produced
#Note there are a number of NAs created due the the log of a negative being not computible##look at lnRR by year
MLMA_2009_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2009)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2009_lnRR)##
## Multivariate Meta-Analysis Model (k = 19; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.2698 90.5396 94.5396 96.3203 95.3396
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 8.8874 2.9812 19 no obs
##
## Test for Heterogeneity:
## Q(df = 18) = 119881336.7140, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 1.7731 0.6867 2.5822 0.0098 0.4273 3.1190 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2010_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2010)
summary(MLMA_2010_lnRR)##
## Multivariate Meta-Analysis Model (k = 16; method: REML)
##
## logLik Deviance AIC BIC AICc
## -34.5830 69.1660 73.1660 74.5821 74.1660
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 5.6352 2.3739 16 no obs
##
## Test for Heterogeneity:
## Q(df = 15) = 64897.7678, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 3.6777 0.6010 6.1197 <.0001 2.4998 4.8556 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2011_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2011)
summary(MLMA_2011_lnRR)##
## Multivariate Meta-Analysis Model (k = 21; method: REML)
##
## logLik Deviance AIC BIC AICc
## -24.6248 49.2495 53.2495 55.2410 53.9554
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1397 0.3738 21 no obs
##
## Test for Heterogeneity:
## Q(df = 20) = 98.8893, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0188 0.1062 0.1771 0.8594 -0.1894 0.2270
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2012_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2012)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2012_lnRR)##
## Multivariate Meta-Analysis Model (k = 47; method: REML)
##
## logLik Deviance AIC BIC AICc
## -91.1326 182.2651 186.2651 189.9224 186.5442
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.8575 1.3629 47 no obs
##
## Test for Heterogeneity:
## Q(df = 46) = 15391.4935, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1019 0.2214 0.4604 0.6452 -0.3320 0.5359
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2013_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2013)
summary(MLMA_2013_lnRR)##
## Multivariate Meta-Analysis Model (k = 54; method: REML)
##
## logLik Deviance AIC BIC AICc
## -104.0133 208.0266 212.0266 215.9671 212.2666
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.5594 1.2487 54 no obs
##
## Test for Heterogeneity:
## Q(df = 53) = 712.3763, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5893 0.2035 -2.8963 0.0038 -0.9880 -0.1905 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2014_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2014)
summary(MLMA_2014_lnRR)##
## Multivariate Meta-Analysis Model (k = 65; method: REML)
##
## logLik Deviance AIC BIC AICc
## -147.9047 295.8093 299.8093 304.1271 300.0060
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 5.8596 2.4207 65 no obs
##
## Test for Heterogeneity:
## Q(df = 64) = 1197354.2596, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2053 0.3031 -0.6773 0.4982 -0.7994 0.3888
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2015_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2015)
summary(MLMA_2015_lnRR)##
## Multivariate Meta-Analysis Model (k = 99; method: REML)
##
## logLik Deviance AIC BIC AICc
## -74.9918 149.9836 153.9836 159.1535 154.1099
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0943 0.3071 99 no obs
##
## Test for Heterogeneity:
## Q(df = 98) = 367.3613, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.1021 0.0429 -2.3811 0.0173 -0.1861 -0.0181 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2016_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2016)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2016_lnRR)##
## Multivariate Meta-Analysis Model (k = 51; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.0134 80.0267 84.0267 87.8508 84.2820
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.2381 0.4880 51 no obs
##
## Test for Heterogeneity:
## Q(df = 50) = 70868.7377, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0421 0.0726 -0.5801 0.5619 -0.1844 0.1002
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2017_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2017)
summary(MLMA_2017_lnRR)##
## Multivariate Meta-Analysis Model (k = 72; method: REML)
##
## logLik Deviance AIC BIC AICc
## -49.9959 99.9917 103.9917 108.5171 104.1682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1058 0.3252 72 no obs
##
## Test for Heterogeneity:
## Q(df = 71) = 255.4688, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0337 0.0480 -0.7033 0.4819 -0.1277 0.0603
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2018_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2018)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2018_lnRR)##
## Multivariate Meta-Analysis Model (k = 148; method: REML)
##
## logLik Deviance AIC BIC AICc
## -245.1928 490.3856 494.3856 500.3665 494.4690
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.9890 0.9945 148 no obs
##
## Test for Heterogeneity:
## Q(df = 147) = 2047.3810, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0073 0.0889 -0.0822 0.9344 -0.1816 0.1670
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MLMA_2019_lnRR <- rma.mv(yi ~ 1, V = vi, random=~1|obs, data=lnRR2019)
summary(MLMA_2019_lnRR)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -47.0733 94.1465 98.1465 102.0102 98.3965
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1813 0.4258 52 no obs
##
## Test for Heterogeneity:
## Q(df = 51) = 302.1940, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0709 0.0646 -1.0963 0.2730 -0.1976 0.0558
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##set prior
prior <- list(R=list(V = 1, nu =0.002), G = list(G = list(V=1, nu = 0.002)))##run bayesian MLMA models
model_magnitude_bayes_2009 <- MCMCglmm(yi ~ 1, mev = lnRR2009clean$vi, random = ~obs, data = lnRR2009clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2010 <- MCMCglmm(yi ~ 1, mev = lnRR2010clean$vi, random = ~obs, data = lnRR2010clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2011 <- MCMCglmm(yi ~ 1, mev = lnRR2011clean$vi, random = ~obs, data = lnRR2011clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2012 <- MCMCglmm(yi ~ 1, mev = lnRR2012clean$vi, random = ~obs, data = lnRR2012clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2013 <- MCMCglmm(yi ~ 1, mev = lnRR2013clean$vi, random = ~obs, data = lnRR2013clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2014 <- MCMCglmm(yi ~ 1, mev = lnRR2014clean$vi, random = ~obs, data = lnRR2014clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2015 <- MCMCglmm(yi ~ 1, mev = lnRR2015clean$vi, random = ~obs, data = lnRR2015clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2016 <- MCMCglmm(yi ~ 1, mev = lnRR2016clean$vi, random = ~obs, data = lnRR2016clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2017 <- MCMCglmm(yi ~ 1, mev = lnRR2017clean$vi, random = ~obs, data = lnRR2017clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2018 <- MCMCglmm(yi ~ 1, mev = lnRR2018clean$vi, random = ~obs, data = lnRR2018clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)
model_magnitude_bayes_2019 <- MCMCglmm(yi ~ 1, mev = lnRR2019clean$vi, random = ~obs, data = lnRR2019clean, prior = prior, burnin = 10000, nitt = 1000000, thin = 100, verbose = FALSE)##get model summaries
summary(model_magnitude_bayes_2009)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 31.51985
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 5.097 0.0002356 14.26 3103
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 5.061 0.0002329 14.29 3114
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 1.7739 0.3589 3.2469 9536 0.0208 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_magnitude_bayes_2010)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 23.90192
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 3.341 0.0001892 9.677 4585
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 3.338 0.0002461 9.757 4192
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.673 2.425 4.973 9900 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_magnitude_bayes_2011)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -11.28967
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.08409 0.0002079 0.2392 8200
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.08199 0.0002033 0.2395 7772
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.01837 -0.21101 0.23836 9900 0.885summary(model_magnitude_bayes_2012)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 45.14502
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.9774 0.0001899 2.489 1733
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1.006 0.0001908 2.499 1713
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.09536 -0.37585 0.53135 9900 0.669summary(model_magnitude_bayes_2013)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 39.93272
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.8623 0.0002047 2.075 1589
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.785 0.0001698 2.025 1654
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.5876 -0.9980 -0.1784 9900 0.00667 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_magnitude_bayes_2014)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 87.08115
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 3.105 0.0001631 7.274 850.3
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 2.965 0.0001571 7.231 834.2
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.2140 -0.8044 0.4100 9900 0.487summary(model_magnitude_bayes_2015)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -68.60624
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.04967 0.0002081 0.1198 2438
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.04949 0.0002497 0.1195 2514
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.10226 -0.19303 -0.01937 9900 0.0212 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_magnitude_bayes_2016)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -16.68832
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.1268 0.0002429 0.3028 3099
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.1239 0.0001545 0.3033 3152
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.04269 -0.18783 0.10043 9900 0.56summary(model_magnitude_bayes_2017)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -47.45125
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.05514 0.0001951 0.1385 3296
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.05638 0.0001962 0.1418 3313
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.03344 -0.13233 0.05938 9900 0.487summary(model_magnitude_bayes_2018)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 88.38417
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.5 0.0002256 1.146 604.6
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.5085 0.0002445 1.153 625.3
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.008286 -0.184830 0.165156 9900 0.928summary(model_magnitude_bayes_2019)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -23.5312
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.09752 0.0002232 0.2482 3376
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.09594 0.0001801 0.2484 3712
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.07209 -0.20128 0.05897 9900 0.269##extract posteriors
sol2009 <- model_magnitude_bayes_2009$Sol
VCV2009 <- model_magnitude_bayes_2009$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2009$VCV))]
sol2010 <- model_magnitude_bayes_2010$Sol
VCV2010 <- model_magnitude_bayes_2010$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2010$VCV))]
sol2011 <- model_magnitude_bayes_2011$Sol
VCV2011 <- model_magnitude_bayes_2011$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2011$VCV))]
sol2012 <- model_magnitude_bayes_2012$Sol
VCV2012 <- model_magnitude_bayes_2012$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2012$VCV))]
sol2013 <- model_magnitude_bayes_2013$Sol
VCV2013 <- model_magnitude_bayes_2013$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2013$VCV))]
sol2014 <- model_magnitude_bayes_2014$Sol
VCV2014 <- model_magnitude_bayes_2014$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2014$VCV))]
sol2015 <- model_magnitude_bayes_2015$Sol
VCV2015 <- model_magnitude_bayes_2015$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2015$VCV))]
sol2016 <- model_magnitude_bayes_2016$Sol
VCV2016 <- model_magnitude_bayes_2016$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2016$VCV))]
sol2017 <- model_magnitude_bayes_2017$Sol
VCV2017 <- model_magnitude_bayes_2017$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2017$VCV))]
sol2018 <- model_magnitude_bayes_2018$Sol
VCV2018 <- model_magnitude_bayes_2018$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2018$VCV))]
sol2019 <- model_magnitude_bayes_2019$Sol
VCV2019 <- model_magnitude_bayes_2019$VCV[,-match("sqrt(mev):sqrt(mev).meta", colnames(model_magnitude_bayes_2019$VCV))]##get folded normal function
mu.fnorm <- function(mu, sigma){dnorm(mu, 0, sigma)*2*sigma^2 + mu*(2*pnorm(mu, 0, sigma) -1)}##get magnitude means + variance
magnitude_mean_2009 <- mu.fnorm(sol2009[,1], sqrt(rowSums(VCV2009)))
magnitude_2009_mean <- data.frame(mean_mag = mean(magnitude_mean_2009), L_CI = HPDinterval(magnitude_mean_2009)[1], U_CI = HPDinterval(magnitude_mean_2009)[2])
magnitude_mean_2010 <- mu.fnorm(sol2010[,1], sqrt(rowSums(VCV2010)))
magnitude_2010_mean <- data.frame(mean_mag = mean(magnitude_mean_2010), L_CI = HPDinterval(magnitude_mean_2010)[1], U_CI = HPDinterval(magnitude_mean_2010)[2])
magnitude_mean_2011 <- mu.fnorm(sol2011[,1], sqrt(rowSums(VCV2011)))
magnitude_2011_mean <- data.frame(mean_mag = mean(magnitude_mean_2011), L_CI = HPDinterval(magnitude_mean_2011)[1], U_CI = HPDinterval(magnitude_mean_2011)[2])
magnitude_mean_2012 <- mu.fnorm(sol2012[,1], sqrt(rowSums(VCV2012)))
magnitude_2012_mean <- data.frame(mean_mag = mean(magnitude_mean_2012), L_CI = HPDinterval(magnitude_mean_2012)[1], U_CI = HPDinterval(magnitude_mean_2012)[2])
magnitude_mean_2013 <- mu.fnorm(sol2013[,1], sqrt(rowSums(VCV2013)))
magnitude_2013_mean <- data.frame(mean_mag = mean(magnitude_mean_2013), L_CI = HPDinterval(magnitude_mean_2013)[1], U_CI = HPDinterval(magnitude_mean_2013)[2])
magnitude_mean_2014 <- mu.fnorm(sol2014[,1], sqrt(rowSums(VCV2014)))
magnitude_2014_mean <- data.frame(mean_mag = mean(magnitude_mean_2014), L_CI = HPDinterval(magnitude_mean_2014)[1], U_CI = HPDinterval(magnitude_mean_2014)[2])
magnitude_mean_2015 <- mu.fnorm(sol2015[,1], sqrt(rowSums(VCV2015)))
magnitude_2015_mean <- data.frame(mean_mag = mean(magnitude_mean_2015), L_CI = HPDinterval(magnitude_mean_2015)[1], U_CI = HPDinterval(magnitude_mean_2015)[2])
magnitude_mean_2016 <- mu.fnorm(sol2016[,1], sqrt(rowSums(VCV2016)))
magnitude_2016_mean <- data.frame(mean_mag = mean(magnitude_mean_2016), L_CI = HPDinterval(magnitude_mean_2016)[1], U_CI = HPDinterval(magnitude_mean_2016)[2])
magnitude_mean_2017 <- mu.fnorm(sol2017[,1], sqrt(rowSums(VCV2017)))
magnitude_2017_mean <- data.frame(mean_mag = mean(magnitude_mean_2017), L_CI = HPDinterval(magnitude_mean_2017)[1], U_CI = HPDinterval(magnitude_mean_2017)[2])
magnitude_mean_2018 <- mu.fnorm(sol2018[,1], sqrt(rowSums(VCV2018)))
magnitude_2018_mean <- data.frame(mean_mag = mean(magnitude_mean_2018), L_CI = HPDinterval(magnitude_mean_2018)[1], U_CI = HPDinterval(magnitude_mean_2018)[2])
magnitude_mean_2019 <- mu.fnorm(sol2019[,1], sqrt(rowSums(VCV2019)))
magnitude_2019_mean <- data.frame(mean_mag = mean(magnitude_mean_2019), L_CI = HPDinterval(magnitude_mean_2019)[1], U_CI = HPDinterval(magnitude_mean_2019)[2])##view ES magnitudes and uncertainty
magnitude_2009_mean## mean_mag L_CI U_CI
## 1 2.958831 2.049309 3.96643magnitude_2010_mean## mean_mag L_CI U_CI
## 1 3.894468 2.84035 4.952511magnitude_2011_mean## mean_mag L_CI U_CI
## 1 0.3286311 0.1887599 0.4890088magnitude_2012_mean## mean_mag L_CI U_CI
## 1 1.132439 0.8757109 1.42672magnitude_2013_mean## mean_mag L_CI U_CI
## 1 1.135282 0.8670486 1.406788magnitude_2014_mean## mean_mag L_CI U_CI
## 1 1.980475 1.641251 2.347749magnitude_2015_mean## mean_mag L_CI U_CI
## 1 0.2648629 0.2009352 0.3320272magnitude_2016_mean## mean_mag L_CI U_CI
## 1 0.4023895 0.3132693 0.5012349magnitude_2017_mean## mean_mag L_CI U_CI
## 1 0.2675296 0.1901038 0.3501495magnitude_2018_mean## mean_mag L_CI U_CI
## 1 0.8021178 0.6840414 0.925188magnitude_2019_mean## mean_mag L_CI U_CI
## 1 0.3553505 0.2536517 0.4683622#Note Code below not in original. Made to allow a data frame to be constructed and plotted from.
magnitudedata = rbind(magnitude_2009_mean, magnitude_2010_mean, magnitude_2011_mean, magnitude_2012_mean, magnitude_2013_mean, magnitude_2014_mean, magnitude_2015_mean, magnitude_2016_mean, magnitude_2017_mean, magnitude_2018_mean, magnitude_2019_mean)yearlabel = c("2009", "2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019")magnitudedata = cbind(yearlabel, magnitudedata )
magnitudedata## yearlabel mean_mag L_CI U_CI
## 1 2009 2.9588308 2.0493088 3.9664297
## 2 2010 3.8944684 2.8403496 4.9525107
## 3 2011 0.3286311 0.1887599 0.4890088
## 4 2012 1.1324393 0.8757109 1.4267201
## 5 2013 1.1352823 0.8670486 1.4067885
## 6 2014 1.9804749 1.6412511 2.3477486
## 7 2015 0.2648629 0.2009352 0.3320272
## 8 2016 0.4023895 0.3132693 0.5012349
## 9 2017 0.2675296 0.1901038 0.3501495
## 10 2018 0.8021178 0.6840414 0.9251880
## 11 2019 0.3553505 0.2536517 0.4683622#Note Plot from above model was not included in the original code. This has been built from scratch. Should be the same as Fig 1B
Decline_magnitude<-ggplot(magnitudedata,aes(x=yearlabel, y=mean_mag, colour=yearlabel)) + geom_line(aes(group=1)) + scale_color_viridis(discrete=TRUE)+ geom_point(size=3) + geom_errorbar(aes (ymin = L_CI, ymax = U_CI), width=0.2) + scale_y_continuous(breaks = seq(0, 15, by = 1), minor_breaks = NULL, limits=c(0,8),labels = scales::number_format(accuracy = 0.1)) + theme(legend.position = "none") + xlab("Year of publication online") + ylab("Effect size magnitude (lnRR)") + theme_minimal(12) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme(legend.position = "none") + theme(axis.ticks=element_line()) + theme(axis.text.x = element_text(angle = 45, hjust=1))+ geom_hline(yintercept = 0)
Decline_magnitude
ggsave(Decline_magnitude, filename = 'Decline magnitude corrected and complete 0.1 and 0.001 remove CO2 under 800 model.png', device=png, width = 4.2, height = 4.6, units = "in", res = 800)CREATE SCATTERPLOT FIGURE TO VISUALIZE MEAN EFFECT SIZE MAGNITUDE FOR EACH OBSERVATION OVER TIME (FIG 1A)
##attach dataset
decline_allobs<-read.csv(file.choose()) ##use dataset "S10 Data"
attach(decline_allobs)## The following objects are masked from decline:
##
## ctrl.n, obs, study, year.online, year.printsummary(decline_allobs)## obs study year.online year.print
## Min. : 1.0 Length:839 Min. :2009 Min. :2009
## 1st Qu.:210.5 Class :character 1st Qu.:2012 1st Qu.:2013
## Median :420.0 Mode :character Median :2015 Median :2015
## Mean :420.0 Mean :2015 Mean :2015
## 3rd Qu.:629.5 3rd Qu.:2017 3rd Qu.:2018
## Max. :839.0 Max. :2019 Max. :2019
##
## ctrl.n lnrr.mag
## Min. : 3.00 Min. :0.0000
## 1st Qu.: 10.00 1st Qu.:0.1021
## Median : 18.00 Median :0.2738
## Mean : 28.77 Mean :0.8639
## 3rd Qu.: 30.00 3rd Qu.:0.8289
## Max. :752.00 Max. :8.4360
## NA's :195##Create plot
#Note I had to fix the specifications of the scale continuous both x and y Original below not working: scale_x_continuous(breaks = round(seq(min(study\(Year), max(study\)Year), by = 1),1)) + scale_y_continuous(breaks = round(seq(min(study$lnrr), 15, by = 1),1)) #Note Code given did not produce a graph that looked exactly the same as in the figures of the paper. Changes to aesthetics were add to make it look the same as the published version.
Decline_studies_loess<-ggplot(decline_allobs,aes(x=year.online, y=lnrr.mag, color=study)) + geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95,color="black") + geom_point(size=ctrl.n*0.03,alpha=0.6) + scale_size(range = c(1, 2), name="Sample size")+ scale_color_viridis(discrete=TRUE)+ xlab("Year of publication online")+ylab("Effect size magnitude (lnRR)") + scale_x_continuous(breaks = round(seq(min(2009), max(2019), by = 1),1)) + scale_y_continuous(breaks = seq(0, 14, by = 1), minor_breaks = NULL, limits=c(-1,14),labels = scales::number_format(accuracy = 0.1)) + theme_minimal(12) + theme(legend.position = "none") + theme(panel.grid.minor = element_blank()) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme(legend.position = "none")+ theme(axis.ticks=element_line()) + theme(axis.text.x = element_text(angle = 45, hjust=1))+ geom_hline(yintercept = 0)
Decline_studies_loess## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 195 rows containing non-finite values (stat_smooth).## Warning: Removed 195 rows containing missing values (geom_point).
ggsave(Decline_studies_loess, filename = 'Decline magnitude corrected and complete 0.1 and 0.001 remove CO2 under 800 lnRR.png', device=png, width = 4.2, height = 4.6, units = "in", res = 800)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 195 rows containing non-finite values (stat_smooth).## Warning: Removed 195 rows containing missing values (geom_point).