Rerun of Clements et al data original
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)#Note I noticed there is an issues in the S5 file published online as the file that was use. There is an issue with year of publishing and year online. This different to the full excel file published online. Without fixing the R code cannot run 2009 as there is nothing in 2009. The S5 year of publishing and online also mismatches with the S10 file for creating Fig 1A. New S5.csv file created from full excel file published
head(decline)## obs study authors year.online year.print if.at.pub X2017.if if.group
## 1 1 a1 Munday et al 2009 2009 9.432 9.504 J
## 2 2 a1 Munday et al 2009 2009 9.432 9.504 J
## 3 3 a1 Munday et al 2009 2009 9.432 9.504 J
## 4 4 a1 Munday et al 2009 2009 9.432 9.504 J
## 5 5 a1 Munday et al 2009 2009 9.432 9.504 J
## 6 6 a1 Munday et al 2009 2009 9.432 9.504 J
## avg.n species climate cue cue.type life.stage ctrl.n ctrl.mean
## 1 27 Amphiprion percula Trop Yes Habitat Larvae 26 94.129
## 2 NA Amphiprion percula Trop Yes Habitat Larvae 20 0.783
## 3 NA Amphiprion percula Trop Yes Habitat Larvae 20 46.380
## 4 NA Amphiprion percula Trop Yes Habitat Larvae 10 98.826
## 5 NA Amphiprion percula Trop Yes Kin Larvae 30 0.912
## 6 NA Amphiprion percula Trop Yes Kin Larvae 30 90.876
## ctrl.sd oa.n oa.mean oa.sd
## 1 2.9931 46 72.789 13.2730
## 2 3.5017 46 83.562 10.6211
## 3 14.0023 46 71.429 13.2730
## 4 1.8563 16 88.258 7.0440
## 5 2.9960 20 99.818 0.0045
## 6 3.9984 20 99.818 0.0045##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:818 2018 :153 2018 :175
## 2 : 1 a87 : 40 Class :character 2015 :108 2015 : 98
## 3 : 1 a11 : 39 Mode :character 2017 :101 2012 : 89
## 4 : 1 a90 : 36 2012 : 97 2017 : 85
## 5 : 1 a22 : 30 2014 : 83 2016 : 84
## 6 : 1 a31 : 28 2013 : 58 2013 : 78
## (Other):812 (Other):597 (Other):218 (Other):209
## if.at.pub X2017.if if.group avg.n
## Length:818 Length:818 Length:818 Min. : 4.00
## Class :character Class :character Class :character 1st Qu.: 12.50
## Mode :character Mode :character Mode :character Median : 18.00
## Mean : 32.62
## 3rd Qu.: 30.00
## Max. :568.00
## NA's :727
## species climate cue cue.type
## Length:818 Length:818 Length:818 Length:818
## 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:818 Min. : 4 Min. : -69.78 Min. : 0.000
## Class :character 1st Qu.: 10 1st Qu.: 1.03 1st Qu.: 1.176
## Mode :character Median : 18 Median : 9.20 Median : 5.560
## Mean : 29 Mean : 438.64 Mean : 104.430
## 3rd Qu.: 30 3rd Qu.: 44.01 3rd Qu.: 21.640
## Max. :752 Max. :154936.88 Max. :25490.446
##
## oa.n oa.mean oa.sd
## Min. : 2.00 Min. : -59.67 Min. : 0.00
## 1st Qu.: 11.00 1st Qu.: 1.14 1st Qu.: 1.15
## Median : 18.00 Median : 12.12 Median : 7.10
## Mean : 28.93 Mean : 438.34 Mean : 114.41
## 3rd Qu.: 30.75 3rd Qu.: 43.51 3rd Qu.: 22.43
## 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)##comupte 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)## Warning in log(m1i/m2i): NaNs producedlnRR2012 <- 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)## Warning in log(m1i/m2i): NaNs producedlnRR2013 <- 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)## Warning in log(m1i/m2i): NaNs producedlnRR2014 <- 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)## Warning in log(m1i/m2i): NaNs producedlnRR2016 <- 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)## Warning in log(m1i/m2i): NaNs producedlnRR2018 <- 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: Rows with NAs omitted from model fitting.summary(MLMA_2011_lnRR)##
## Multivariate Meta-Analysis Model (k = 48; method: REML)
##
## logLik Deviance AIC BIC AICc
## -37.4820 74.9639 78.9639 82.6642 79.2367
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1496 0.3868 48 no obs
##
## Test for Heterogeneity:
## Q(df = 47) = 228.8028, 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: Rows with NAs omitted from model fitting.
## 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 = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -149.7127 299.4253 303.4253 308.3572 303.5682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.2064 1.0984 88 no obs
##
## Test for Heterogeneity:
## Q(df = 87) = 29290.8615, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1353 0.1308 1.0343 0.3010 -0.1211 0.3917
##
## ---
## 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: Rows with NAs omitted from model fitting.
## 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 = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -102.2558 204.5116 208.5116 212.3752 208.7616
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.2182 1.1037 52 no obs
##
## Test for Heterogeneity:
## Q(df = 51) = 321.2026, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5545 0.1848 -3.0002 0.0027 -0.9168 -0.1923 **
##
## ---
## 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)## Warning: Rows with NAs omitted from model fitting.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
## -44.9446 89.8892 93.8892 97.9399 94.1156
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.2268 0.4762 57 no obs
##
## Test for Heterogeneity:
## Q(df = 56) = 70885.9993, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0193 0.0687 -0.2806 0.7790 -0.1539 0.1154
##
## ---
## 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)## Warning: Rows with NAs omitted from model fitting.summary(MLMA_2017_lnRR)##
## Multivariate Meta-Analysis Model (k = 99; method: REML)
##
## logLik Deviance AIC BIC AICc
## -74.2817 148.5633 152.5633 157.7333 152.6897
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1508 0.3884 99 no obs
##
## Test for Heterogeneity:
## Q(df = 98) = 2575.9384, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0320 0.0460 0.6952 0.4869 -0.0582 0.1221
##
## ---
## 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
## -257.9394 515.8788 519.8788 525.9266 519.9594
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.2818 0.5309 153 no obs
##
## Test for Heterogeneity:
## Q(df = 152) = 3035.7779, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0129 0.0493 0.2613 0.7939 -0.0838 0.1095
##
## ---
## 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 = 43; method: REML)
##
## logLik Deviance AIC BIC AICc
## 18.8145 -37.6289 -33.6289 -30.1536 -33.3212
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0097 0.0986 43 no obs
##
## Test for Heterogeneity:
## Q(df = 42) = 100.4464, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0207 0.0207 -1.0004 0.3171 -0.0612 0.0198
##
## ---
## 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: 44.13777
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 17.44 0.0001568 48.92 2364
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 17.43 0.0001937 49.56 2373
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.1533 0.4633 5.8532 9503 0.0222 *
## ---
## 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: 106.4366
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 13.76 0.0002376 33.47 887.5
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 13.52 0.0002257 33.66 914.4
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.699 2.207 5.173 10465 <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: -19.92756
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.08014 0.0001852 0.2064 4213
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.08163 0.0001539 0.208 4372
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.03973 -0.18120 0.10720 9900 0.569summary(model_magnitude_bayes_2012)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 70.82581
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.6042 0.0001614 1.451 1051
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.642 0.0002234 1.463 1050
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.1339 -0.1298 0.3998 9900 0.315summary(model_magnitude_bayes_2013)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: 38.47168
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.6507 0.0002146 1.725 2166
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.6469 0.0002343 1.719 2034
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.5585 -0.9506 -0.1983 9900 0.00263 **
## ---
## 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: 144.7408
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 7.377 0.0002255 16.78 529.8
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 6.819 0.0002211 16.59 537.8
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.08116 -0.91233 0.71634 9900 0.852summary(model_magnitude_bayes_2015)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -72.13293
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.0466 0.0001863 0.1115 2519
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.04671 0.0002183 0.1114 2430
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.09319 -0.17720 -0.01339 9900 0.0269 *
## ---
## 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: -10.84778
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.115 0.0001644 0.2816 2965
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.1234 0.0001433 0.289 2910
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.0195 -0.1579 0.1180 9102 0.779summary(model_magnitude_bayes_2017)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -57.17552
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.08063 0.0003114 0.1843 2012
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.07592 0.0001879 0.1814 1814
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.03215 -0.06386 0.12241 9900 0.497summary(model_magnitude_bayes_2018)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -6.007702
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.1535 0.0002461 0.3968 1150
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.158 0.000209 0.4003 1157
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 0.01318 -0.09012 0.11291 9900 0.796summary(model_magnitude_bayes_2019)##
## Iterations = 10001:999901
## Thinning interval = 100
## Sample size = 9900
##
## DIC: -96.79821
##
## G-structure: ~obs
##
## post.mean l-95% CI u-95% CI eff.samp
## obs 0.005841 0.000214 0.01462 9900
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.005866 0.0002217 0.01489 10170
##
## Location effects: yi ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -0.02155 -0.06700 0.02144 10420 0.329##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.420892 3.831258 7.45244magnitude_2010_mean## mean_mag L_CI U_CI
## 1 5.195638 4.184835 6.219871magnitude_2011_mean## mean_mag L_CI U_CI
## 1 0.3236697 0.2267376 0.4298651magnitude_2012_mean## mean_mag L_CI U_CI
## 1 0.9000583 0.7517323 1.063419magnitude_2013_mean## mean_mag L_CI U_CI
## 1 1.017335 0.7151714 1.335721magnitude_2014_mean## mean_mag L_CI U_CI
## 1 3.015001 2.525651 3.50166magnitude_2015_mean## mean_mag L_CI U_CI
## 1 0.2554101 0.194948 0.3189439magnitude_2016_mean## mean_mag L_CI U_CI
## 1 0.3912036 0.3018465 0.481828magnitude_2017_mean## mean_mag L_CI U_CI
## 1 0.3171702 0.2518118 0.3831636magnitude_2018_mean## mean_mag L_CI U_CI
## 1 0.442555 0.3266386 0.5726379magnitude_2019_mean## mean_mag L_CI U_CI
## 1 0.08786796 0.05404056 0.1281328#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.42089165 3.83125798 7.4524396
## 2 2010 5.19563844 4.18483549 6.2198706
## 3 2011 0.32366974 0.22673764 0.4298651
## 4 2012 0.90005826 0.75173234 1.0634189
## 5 2013 1.01733537 0.71517137 1.3357205
## 6 2014 3.01500136 2.52565066 3.5016596
## 7 2015 0.25541015 0.19494802 0.3189439
## 8 2016 0.39120361 0.30184653 0.4818280
## 9 2017 0.31717016 0.25181175 0.3831636
## 10 2018 0.44255497 0.32663856 0.5726379
## 11 2019 0.08786796 0.05404056 0.1281328#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 original model lnRR.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, study, year.online, year.printsummary(decline_allobs)## ï..obs study year.print year.online
## Min. : 1.0 Length:795 Min. :2009 Min. :2009
## 1st Qu.:211.5 Class :character 1st Qu.:2013 1st Qu.:2012
## Median :419.0 Mode :character Median :2015 Median :2015
## Mean :414.8 Mean :2015 Mean :2015
## 3rd Qu.:619.5 3rd Qu.:2018 3rd Qu.:2017
## Max. :818.0 Max. :2019 Max. :2019
## ctrl.n lnrr.mag
## Min. : 4.00 Min. : 0.00000
## 1st Qu.: 10.00 1st Qu.: 0.09394
## Median : 18.00 Median : 0.26889
## Mean : 29.36 Mean : 0.96125
## 3rd Qu.: 30.00 3rd Qu.: 0.75979
## Max. :752.00 Max. :13.81551##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'
ggsave(Decline_studies_loess, filename = 'Decline magnitude original lnRR.png', device=png, width = 4.2, height = 4.6, units = "in", res = 800)## `geom_smooth()` using formula 'y ~ x'