Rerun of Clements et al data corrected and complete
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 : 53 Length:839 2018 :153 2018 :173
## 2 : 1 a3 : 48 Class :character 2015 :105 2015 : 95
## 3 : 1 a87 : 40 Mode :character 2017 : 85 2016 : 82
## 4 : 1 a90 : 36 2014 : 83 2013 : 78
## 5 : 1 a31 : 28 2012 : 79 2012 : 71
## (Other):781 a22 : 24 (Other):281 (Other):287
## NA's : 53 (Other):610 NA's : 53 NA's : 53
## if.at.pub X2017.if if.group avg.n
## Length:839 Length:839 Length:839 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 :744
## species climate cue cue.type
## Length:839 Length:839 Length:839 Length:839
## 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:839 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
## NA's :53 NA's :53 NA's :53
## 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
## NA's :53 NA's :53 NA's :53##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
## -56.5722 113.1444 117.1444 118.9251 117.9444
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 30.6312 5.5345 19 no obs
##
## Test for Heterogeneity:
## Q(df = 18) = 10724056.2936, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 3.1705 1.2759 2.4849 0.0130 0.6698 5.6713 *
##
## ---
## 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 = 48; method: REML)
##
## logLik Deviance AIC BIC AICc
## -143.4474 286.8948 290.8948 294.5951 291.1675
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 25.9986 5.0989 48 no obs
##
## Test for Heterogeneity:
## Q(df = 47) = 20986.5067, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 3.7009 0.7397 5.0033 <.0001 2.2511 5.1507 ***
##
## ---
## 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)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2011_lnRR)##
## Multivariate Meta-Analysis Model (k = 51; method: REML)
##
## logLik Deviance AIC BIC AICc
## -58.8621 117.7241 121.7241 125.5482 121.9794
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1496 0.3868 51 no obs
##
## Test for Heterogeneity:
## Q(df = 50) = 228.8036, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0398 0.0716 -0.5553 0.5787 -0.1801 0.1006
##
## ---
## 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
## -154.4766 308.9533 312.9533 317.6667 313.1133
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.2983 1.1394 79 no obs
##
## Test for Heterogeneity:
## Q(df = 78) = 29286.2472, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1296 0.1402 0.9247 0.3551 -0.1451 0.4043
##
## ---
## 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)## Warning: Ratio of largest to smallest sampling variance extremely large. May not
## be able to obtain stable results.summary(MLMA_2013_lnRR)##
## Multivariate Meta-Analysis Model (k = 54; method: REML)
##
## logLik Deviance AIC BIC AICc
## -140.2671 280.5342 284.5342 288.4748 284.7742
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.2604 1.1227 54 no obs
##
## Test for Heterogeneity:
## Q(df = 53) = 357.9118, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5433 0.1865 -2.9138 0.0036 -0.9088 -0.1779 **
##
## ---
## 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
## -226.5561 453.1122 457.1122 461.9256 457.2641
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 13.8435 3.7207 83 no obs
##
## Test for Heterogeneity:
## Q(df = 82) = 6974.5766, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0765 0.4110 -0.1860 0.8524 -0.8821 0.7292
##
## ---
## 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.7140 155.4281 159.4281 164.7169 159.5469
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0887 0.2979 105 no obs
##
## Test for Heterogeneity:
## Q(df = 104) = 368.7121, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0936 0.0410 -2.2807 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) = 70914.5911, 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.0473 122.0946 126.0946 130.9562 126.2427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1465 0.3828 85 no obs
##
## Test for Heterogeneity:
## Q(df = 84) = 2256.3328, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0186 0.0493 0.3776 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
## -331.5337 663.0674 667.0674 673.1152 667.1480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 2.1052 1.4509 153 no obs
##
## Test for Heterogeneity:
## Q(df = 152) = 1394.0520, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0544 0.1245 0.4372 0.6619 -0.1895 0.2983
##
## ---
## 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.1947, 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: 45.55427
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 17.18 0.0002394 49.52 2214
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 17.69 0.0002183 49.04 2376
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.1712 0.5494 5.9057 9900 0.0236 *
## ---
## 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: 104.7543
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 13.64 0.0001605 33.53 891.8
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 13.57 0.0001642 33.82 875.5
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.698 2.291 5.278 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: -23.08193
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.08101 0.0002027 0.2054 3997
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.08043 0.000245 0.2094 4057
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.03939 -0.18108 0.11146 9900 0.596summary(model_magnitude_bayes_2012)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 57.65649
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.6677 0.0001364 1.584 1035
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.6745 0.0002474 1.588 1012
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.1285 -0.1510 0.3994 9518 0.356summary(model_magnitude_bayes_2013)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 39.48676
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.6747 0.0002627 1.785 1965
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.6718 0.000167 1.832 1904
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.5414 -0.9413 -0.1791 9900 0.00303 **
## ---
## 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: 151.9536
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 7.253 0.0002048 16.82 537.2
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 6.988 0.0001676 16.83 554.8
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.07792 -0.84119 0.77089 9587 0.851summary(model_magnitude_bayes_2015)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -87.89382
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.04845 0.0001846 0.1135 2587
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.04497 0.0002298 0.1113 2465
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.093281 -0.172364 -0.009261 10199 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: -16.03626
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.1177 0.0001935 0.2836 2847
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.1186 0.0002459 0.285 2936
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.03658 -0.17055 0.09542 9071 0.587summary(model_magnitude_bayes_2017)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -46.1376
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.07762 0.0002488 0.1792 2313
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.07442 0.0002156 0.1766 2275
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.01845 -0.08129 0.11552 9900 0.71summary(model_magnitude_bayes_2018)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 145.8749
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 1.058 0.0001457 2.575 545.9
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1.072 0.000226 2.563 510.4
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.0551 -0.1899 0.2963 10093 0.659summary(model_magnitude_bayes_2019)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -22.54494
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.09601 0.0002069 0.2497 3544
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.09864 0.0002417 0.251 3490
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.07175 -0.20470 0.05802 9974 0.275##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 5.430454 3.737425 7.317699magnitude_2010_mean## mean_mag L_CI U_CI
## 1 5.19123 4.210107 6.248693magnitude_2011_mean## mean_mag L_CI U_CI
## 1 0.3237188 0.221841 0.4222023magnitude_2012_mean## mean_mag L_CI U_CI
## 1 0.9328352 0.761204 1.101861magnitude_2013_mean## mean_mag L_CI U_CI
## 1 1.026047 0.7209365 1.350764magnitude_2014_mean## mean_mag L_CI U_CI
## 1 3.018606 2.544318 3.540284magnitude_2015_mean## mean_mag L_CI U_CI
## 1 0.2555018 0.194544 0.3185084magnitude_2016_mean## mean_mag L_CI U_CI
## 1 0.3901898 0.3034971 0.4812106magnitude_2017_mean## mean_mag L_CI U_CI
## 1 0.3120894 0.2473594 0.3807806magnitude_2018_mean## mean_mag L_CI U_CI
## 1 1.162974 0.9161666 1.389701magnitude_2019_mean## mean_mag L_CI U_CI
## 1 0.3563277 0.2549305 0.4714779#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 5.4304536 3.7374252 7.3176988
## 2 2010 5.1912305 4.2101067 6.2486934
## 3 2011 0.3237188 0.2218410 0.4222023
## 4 2012 0.9328352 0.7612040 1.1018614
## 5 2013 1.0260469 0.7209365 1.3507637
## 6 2014 3.0186064 2.5443179 3.5402842
## 7 2015 0.2555018 0.1945440 0.3185084
## 8 2016 0.3901898 0.3034971 0.4812106
## 9 2017 0.3120894 0.2473594 0.3807806
## 10 2018 1.1629736 0.9161666 1.3897014
## 11 2019 0.3563277 0.2549305 0.4714779#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 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 : 1.1539
## 3rd Qu.: 30.00 3rd Qu.: 0.8379
## Max. :752.00 Max. :13.8155
## 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 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).