Rerun of Clements et al data corrected and complete 0.1 and 0.001
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 a3 : 48 Length:786 2018 :153 2018 :173
## 2 : 1 a87 : 40 Class :character 2015 :105 2015 : 95
## 3 : 1 a90 : 36 Mode :character 2017 : 85 2016 : 82
## 4 : 1 a31 : 28 2014 : 83 2013 : 78
## 5 : 1 a22 : 24 2012 : 79 2012 : 71
## 6 : 1 a73 : 22 2016 : 57 2017 : 71
## (Other):780 (Other):588 (Other):224 (Other):216
## if.at.pub X2017.if if.group avg.n
## Length:786 Length:786 Length:786 Min. : 4.0
## Class :character Class :character Class :character 1st Qu.: 12.0
## Mode :character Mode :character Mode :character Median : 18.0
## Mean : 32.2
## 3rd Qu.: 30.0
## Max. :568.0
## NA's :691
## species climate cue cue.type
## Length:786 Length:786 Length:786 Length:786
## 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:786 Min. : 3.00 Min. : -69.78 Min. : 0.000
## Class :character 1st Qu.: 10.00 1st Qu.: 1.18 1st Qu.: 1.130
## Mode :character Median : 18.00 Median : 9.98 Median : 5.402
## Mean : 29.12 Mean : 454.55 Mean : 107.674
## 3rd Qu.: 30.00 3rd Qu.: 45.29 3rd Qu.: 21.646
## Max. :752.00 Max. :154936.88 Max. :25490.446
##
## oa.n oa.mean oa.sd
## Min. : 2.00 Min. : -59.67 Min. : 0.00
## 1st Qu.: 10.00 1st Qu.: 1.38 1st Qu.: 1.08
## Median : 18.00 Median : 13.58 Median : 7.02
## Mean : 29.02 Mean : 454.14 Mean : 117.90
## 3rd Qu.: 34.00 3rd Qu.: 44.43 3rd Qu.: 22.00
## Max. :755.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)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2010_lnRR)##
## Multivariate Meta-Analysis Model (k = 48; method: REML)
##
## logLik Deviance AIC BIC AICc
## -108.2735 216.5470 220.5470 224.2473 220.8197
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 5.7551 2.3990 48 no obs
##
## Test for Heterogeneity:
## Q(df = 47) = 3173441.6049, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 2.1335 0.3512 6.0752 <.0001 1.4452 2.8218 ***
##
## ---
## 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 = 51; method: REML)
##
## logLik Deviance AIC BIC AICc
## -51.9766 103.9531 107.9531 111.7771 108.2084
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1496 0.3868 51 no obs
##
## Test for Heterogeneity:
## Q(df = 50) = 228.8395, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0398 0.0716 -0.5558 0.5783 -0.1801 0.1005
##
## ---
## 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 = 79; method: REML)
##
## logLik Deviance AIC BIC AICc
## -136.7398 273.4797 277.4797 282.1931 277.6397
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.3078 1.1436 79 no obs
##
## Test for Heterogeneity:
## Q(df = 78) = 29295.0394, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0982 0.1398 0.7029 0.4821 -0.1757 0.3722
##
## ---
## 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 = 83; method: REML)
##
## logLik Deviance AIC BIC AICc
## -181.8703 363.7406 367.7406 372.5541 367.8925
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 4.8347 2.1988 83 no obs
##
## Test for Heterogeneity:
## Q(df = 82) = 1333008.5505, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2004 0.2441 -0.8209 0.4117 -0.6789 0.2781
##
## ---
## 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 = 105; method: REML)
##
## logLik Deviance AIC BIC AICc
## -77.7151 155.4302 159.4302 164.7190 159.5490
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0887 0.2979 105 no obs
##
## Test for Heterogeneity:
## Q(df = 104) = 368.7129, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0936 0.0410 -2.2806 0.0226 -0.1740 -0.0132 *
##
## ---
## 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 = 57; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.5933 87.1866 91.1866 95.2373 91.4130
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.2241 0.4734 57 no obs
##
## Test for Heterogeneity:
## Q(df = 56) = 70900.5484, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0350 0.0672 -0.5211 0.6023 -0.1666 0.0966
##
## ---
## 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 = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -61.0486 122.0972 126.0972 130.9589 126.2454
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1465 0.3828 85 no obs
##
## Test for Heterogeneity:
## Q(df = 84) = 2256.3271, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0186 0.0493 0.3777 0.7057 -0.0780 0.1152
##
## ---
## 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 = 153; method: REML)
##
## logLik Deviance AIC BIC AICc
## -250.1242 500.2485 504.2485 510.2962 504.3290
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.9423 0.9707 153 no obs
##
## Test for Heterogeneity:
## Q(df = 152) = 2102.9211, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0015 0.0853 -0.0176 0.9859 -0.1688 0.1658
##
## ---
## 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: 33.87178
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 5.027 0.000176 14.28 3167
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 5.109 0.0002194 14.32 3089
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 1.7777 0.3983 3.2976 10705 0.0196 *
## ---
## 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: 65.64687
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 3.063 0.0002331 7.462 1154
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 2.988 0.0001512 7.382 1161
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 2.129 1.415 2.814 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: -20.12467
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.07962 0.0002012 0.2075 3703
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.08182 0.0001623 0.2093 3815
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.04042 -0.18871 0.10239 10422 0.574summary(model_magnitude_bayes_2012)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 52.4577
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.6926 0.0002066 1.616 1121
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.6628 0.0001806 1.591 1107
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.09584 -0.17627 0.38119 10888 0.493summary(model_magnitude_bayes_2013)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 47.76008
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.8277 0.0002142 2.068 1491
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.8246 0.0001592 2.065 1571
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.5950 -1.0012 -0.1878 9900 0.00404 **
## ---
## 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: 118.141
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 2.45 0.0002257 5.788 726.3
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 2.532 0.0001376 5.788 731.6
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.2020 -0.6803 0.2828 9900 0.407summary(model_magnitude_bayes_2015)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -71.07368
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.04608 0.0001961 0.1103 2653
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.04648 0.0002746 0.1099 2538
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.09294 -0.17684 -0.01275 10181 0.0253 *
## ---
## 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: -12.35051
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.1138 0.0002643 0.2813 3012
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.1214 0.0001925 0.2805 2955
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.03498 -0.17145 0.10123 10457 0.601summary(model_magnitude_bayes_2017)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -44.89313
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.07757 0.0001748 0.1787 2171
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.07523 0.0002385 0.1789 2220
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.01788 -0.07618 0.11993 9900 0.711summary(model_magnitude_bayes_2018)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 106.0556
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.4587 0.0001917 1.076 616.2
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.5009 0.0002283 1.094 620.1
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.001203 -0.170626 0.170371 9900 0.983summary(model_magnitude_bayes_2019)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -21.36962
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.0954 0.0002187 0.2468 3537
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.09737 0.0003016 0.2473 3596
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.07170 -0.19618 0.06515 10116 0.27##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.957034 2.059286 3.982292magnitude_2010_mean## mean_mag L_CI U_CI
## 1 2.667045 2.147438 3.162456magnitude_2011_mean## mean_mag L_CI U_CI
## 1 0.3235893 0.2258914 0.4261405magnitude_2012_mean## mean_mag L_CI U_CI
## 1 0.9347203 0.7719635 1.119563magnitude_2013_mean## mean_mag L_CI U_CI
## 1 1.139086 0.8731541 1.428216magnitude_2014_mean## mean_mag L_CI U_CI
## 1 1.79311 1.517364 2.083125magnitude_2015_mean## mean_mag L_CI U_CI
## 1 0.2543803 0.1912708 0.3144037magnitude_2016_mean## mean_mag L_CI U_CI
## 1 0.3892753 0.3026903 0.4784667magnitude_2017_mean## mean_mag L_CI U_CI
## 1 0.3128215 0.2478361 0.3829406magnitude_2018_mean## mean_mag L_CI U_CI
## 1 0.7823563 0.6708826 0.9046148magnitude_2019_mean## mean_mag L_CI U_CI
## 1 0.3546911 0.2520477 0.466605#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.9570336 2.0592856 3.9822925
## 2 2010 2.6670445 2.1474376 3.1624555
## 3 2011 0.3235893 0.2258914 0.4261405
## 4 2012 0.9347203 0.7719635 1.1195629
## 5 2013 1.1390860 0.8731541 1.4282156
## 6 2014 1.7931096 1.5173638 2.0831250
## 7 2015 0.2543803 0.1912708 0.3144037
## 8 2016 0.3892753 0.3026903 0.4784667
## 9 2017 0.3128215 0.2478361 0.3829406
## 10 2018 0.7823563 0.6708826 0.9046148
## 11 2019 0.3546911 0.2520477 0.4666050#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 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.1020
## Median : 18.00 Median :0.2806
## Mean : 28.77 Mean :0.8434
## 3rd Qu.: 30.00 3rd Qu.:0.8379
## Max. :752.00 Max. :8.4360
## NA's :53##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 53 rows containing non-finite values (stat_smooth).## Warning: Removed 53 rows containing missing values (geom_point).
ggsave(Decline_studies_loess, filename = 'Decline magnitude corrected and complete 0.1 and 0.001 lnRR.png', device=png, width = 4.2, height = 4.6, units = "in", res = 800)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 53 rows containing non-finite values (stat_smooth).## Warning: Removed 53 rows containing missing values (geom_point).