Merge branch 'master' of github.com:daniel1noble/ecology_comment

.gitignore

@@ -3,7 +3,6 @@ Ecology\ comments/
 Ecology_Comment/
 DataS1_ORIGINAL_DO_NOT_TOUCH/
 Supplemental_Implementation_Example_cache/
-Supplemental_Implementation_Example.html
 sim_check_rerun/
 *.docx
 New_sims.pdf

README.md

@@ -1 +1,14 @@
-# ecology_comment
+# Ecology Commentary 
+
+## Introduction
+This repository houses the code and supplementary tutorial used to demonstrate how multi-level meta-analytic models from `metafor` can be corrected to avoid infated Type I error in the presence of non-independent effect sizes. The commentary is a response to Song et al. (2020), to show how a few simple corrections can provide some resolution to problems they identify in their very thorough simulations.
+
+*Reproducing the simulations?* Users wanting to reproduce Song et al's (2020) simulations, along with the correlations implemented by us can find all the R code in tge `R/` directory.
+
+*Just want to know how to apply corrections?* Users who are interested in learning more about how they can correct for non-independence can open and read the `Supplemental_Implementation_Example.html`
+
+## Citation
+Nakagawa, S., A. M. Senior, W. Viechtbauer, and D. W. A. Noble. 2021. An assessment of statistical methods for non-independent data in ecological meta-analyses: comment. Ecology
+
+## References
+Song, C., S. D. Peacor, C. W. Osenberg, and J. R. Bence. 2020. An assessment of statistical methods for nonindependent data in ecological meta-analyses. Ecology online: e03184.

Supplemental_Implementation_Example.Rmd

@@ -1,5 +1,5 @@
 ---
-  title: "Implementation Examples of multi-level meta-analytic models with robust variance estimators and using Bayesian approaches to deal with non-independent effect sizes"
+  title: "Tutorial on how to implement simple corrections to deal with non-independent effect sizes in multi-level meta-analysis"
   author: Shinichi Nakagawa, Alistair M. Senior, Wolfgang Viechtbauer and Daniel W. A. Noble
   date: "`r Sys.Date()`"
   bibliography: refs.bib
@@ -11,25 +11,20 @@
       toc: yes
       toc_depth: 6
       toc_float: yes
-    bookdown::word_document2:
-      toc: no
-      toc_depth: 6
-      number_sections: true
-      reference_docx: template.docx
-  editor_options: 
-    chunk_output_type: console
 ---
 
 ```{r setup, include=FALSE}
   knitr::opts_chunk$set(echo = TRUE, cache = FALSE, message = FALSE, warning = FALSE, tidy = TRUE, fig.width = 10)
   ## numbers >= 10^5 will be denoted in scientific notation,      ## numbers >= 10^5 will be denoted in scientific notation,
     ## and rounded to 2 digits      ## and rounded to 2 digits
     options(digits = 2)
+    
+    #install.packages("devtools")
+#remotes::install_github("rlesur/klippy")
+#library(klippy)
 ```
 
-```{r klippy, echo=FALSE, include=TRUE,results='hide'}
-install.packages("devtools")
-remotes::install_github("rlesur/klippy")
+```{r klippy, echo=FALSE, include=TRUE}
 klippy::klippy(tooltip_message = 'Click to Copy Code', tooltip_success = 'Done', position = 'right', color = "red")
 ```
 
@@ -39,35 +34,33 @@ In this tutorial, we demonstrate how meta-analyst's can implement approaches for
 # R Packages Required
 
 First, we'll load some of the packages that we'll need.
-
-```{r packages_data,results='hide'}
+```{r packages_data,results='hide', class.source='klippy'}
   # Clean workspace
       rm(list=ls())
 
   # Loading packages & Functions
-  
-      install.packages("pacman")
-      pacman::p_load(tidyverse, MASS, kableExtra, gridExtra,  MCMCglmm, brms,  metafor, robumeta, clubSandwich, pander, tidyverse)
+
+  pacman::p_load(tidyverse, MASS, kableExtra, gridExtra,  MCMCglmm, brms,  metafor, robumeta, clubSandwich, pander, tidyverse)
 ```
 
 # Simulating Non-independent Effect Size Data
 
 Here, we will simulate some meta-analytic data. We will keep this very simple, just for demonstration purposes. Hence, we will assume that we have collected data from a total of 20 studies, and we'll assume that we were able to extract n = 3 effect sizes from each of these 20 studies. In total, we have a data set that contains n = `r 20*3` effect sizes.
 
 ```{r simulated data, class.source='klippy'}
-    # Simulate a dataset composed of 30 papers, each having 5 effects from the same study.
+    # Simulate a dataset composed of 20 papers, each having 3 effects from the same study.
   
         set.seed(87)
 
     # Parameters
           no.paper = 20   # Numbers of unique papers
          n.effects = 3    # Number of effects per paper. We will keep simple and balanced
-             rho.e = 0.8  # Correlation in sampling variances
+             rho.e = 0.8  # Correlation among sampling variances
           study_id = rep(1:no.paper, each = n.effects)  # Study ID's
          var_paper = 1   # Between-study variance (i.e., tau2)
         var_effect = 0.8  # Effect size (or within study) variance
                 mu = 0.4  # Average, or overall, effect size
-               rho = 0.1 #c(0.1, 0.2, 0.3, 0.4) # Note that, setting this to zero equates to Shinichi's approach
+               rho = 0.1  # Correlation among effect sizes within study; could vary
    
      # Add sampling variance
           # First, sample
@@ -99,7 +92,7 @@ Here, we will simulate some meta-analytic data. We will keep this very simple, j
       # Build the full correlation matrix
       cor_yi <- as.matrix(Matrix::bdiag(matrices))
       
-      # Calculate the full covariance matrix with Tau2
+      # Calculate the full covariance matrix 
       cov_yi <- cor_yi * sqrt(var_effect) * sqrt(var_effect)
   
     # Now simulate effect sizes, assuming that the average effect size and all the relevant within study correlations. 
@@ -114,9 +107,9 @@ Here, we will simulate some meta-analytic data. We will keep this very simple, j
    
 ```
 
-Now that we have our simulated (i.e., fake or toy) data, we can demonstrate a few corrections that can be applied to MLMA models that will offset any possible inflated Type I error rates.
+Now that we have our simulated data, we can demonstrate a few corrections that can be applied to MLMA models that will offset any possible inflated Type I error rates.
 
-# Fit the Multi-level Meta-analytic (MLMA) Model
+# Step 1: Fit the Multi-level Meta-analytic (MLMA) Model
 
 First, lets just fit our multilevel meta-analytic (MLMA) model. We can do that using our simulated data as follows:
 
@@ -126,40 +119,32 @@ First, lets just fit our multilevel meta-analytic (MLMA) model. We can do that u
                                     data=data, test="t")
   summary(mod_multilevel)
 ```
-We have fit a simple MLMA model that estimates the overall meta-analytic mean. We can see that our model estimates this correctly. Remember that the true mean is `r mu`, and we are pretty close to this value (i.e., `r mod_multilevel$beta`). This is also true of our random effect variance estimates. In this case, we know that the MLMA model is ignoring the shared sampling variance dependence in effect sizes. We expect that this 'should' (at least on average) inflate Type I error rates. Assuming we did not know any better we would want to account for this dependence. Below, we describe a few corrections that meta-analysts can apply that should overcome the problems associated with not accounting for effect size dependence. 
+We have fit a simple MLMA model that estimates the overall meta-analytic mean. We can see that our model estimates this correctly. Remember that the true mean is `r mu`, and we are pretty close to this value (i.e., `r mod_multilevel$beta`). This is also true of our random effect variance estimates. In this case, we know that the MLMA model is not completely dealing with the dependence among effect sizes. We expect that this 'should' (at least on average) inflate Type I error rates. Assuming we did not know any better we would want to account for this dependence. Below, we describe a few corrections that meta-analysts can apply that should overcome the problems associated with not accounting for effect size dependence. 
 
-# Correction using Robust Variance Estimator (RVE) with Saitterwaite Degrees of Freedom Correction
+# Correction 1: Using a Robust Variance Estimator (RVE) with Saitterwaite Degrees of Freedom Correction
 
-A very simple solution is to make use of robust variance estimators (RVE). This can be done in a few packages, but very easily using the `clubSandwich` package [@Pustejovsky2020; @Hedges2010; @Tipon2015], which also works well with `metafor` [@Wolfgang2010] models. This appaoch also makes use of a Saitterwaite degrees of freedom correction [@SW; @Tipon2015]. This works with `metafor` objects quite elegantly. The benefit of such an approach is simply that we need not make any assumptions about what the correlation between effect sizes actually is (assuming we didn't know the true correlation) [@Hedges2010; @Tipon2015]. In addition, it also will acount for possible heteroscedascity. This solution can be implemented as follows using our MLMA model we fit in the above section.
+A very simple solution is to make use of robust variance estimators (RVE). This can be done in a few packages, but very easily using the `clubSandwich` package [@Pustejovsky2020; @Hedges2010; @Tipon2015], which also works well with `metafor` [@Wolfgang2010] models. This approach also makes use of a Saitterwaite degrees of freedom correction [@SW; @Tipon2015]. This works with `metafor` objects quite elegantly. The benefit of such an approach is simply that we need not make any assumptions about what the correlation between effect sizes actually is (assuming we didn't know the true correlation) [@Hedges2010; @Tipon2015]. In addition, it also will account for possible heteroscedascity. This solution can be implemented as follows using our MLMA model we fit in the above section.
 
 ```{r fitrobust, class.source='klippy'}
    mod_RVE <- coef_test(mod_multilevel, vcov="CR2", cluster = data$study_id)
    print(mod_RVE)
-   
-   mod <- robumeta::robu(formula=yi~1, data=data, studynum=study_id, var.eff.size=vi, method = "HIER", small = FALSE)
+  
 ```
 
-With this simple (and well balanced) data, our RVE approches doesn't change the results much. 
+A better, but slightly more restricted RVE can be implemented in the `robumeta` package in R. It is better at dealing with non-independence, but is currently limited to a single random effect level. Nonetheless, with our simple model we can fit a RVE model that completely deals with non-independence as follows:
 
-# Correction by Applying Bayesian Multi-level Meta-analytic Model
-
-As we describe in our comment, Bayesian approaches, assuming one has a good sample size, do a very good job correcting for inflated Type I errors across a variety of situations. Baysian MLMA models can be fit in various packages. Probably the most flexible for meta-analyst's are `MCMCglmm` [@Hadfield2010] and `brms` [@Brkner2017; @Brkner2018]. Here, we demonstrate how to fit the same MLMA model using `MCMCglmm` which has a syntax that is different from the typical one used in packages such as `metafor` and `lme4` [@Bates2015], which meta-analysts might be more acustomed too. 
-
-```{r bayes, class.source='klippy'}
-    prior <- list(R = list(V = 1, nu = 0.002), 
-                  G = list(G1 = list(V = 1 , nu = 1, alpha.mu=0, alpha.V=25^2)))
-
-   bayes_multilevel <- MCMCglmm(yi ~ 1, mev = data$vi, random = ~ study_id, data = data, prior = prior, verbose = FALSE)
-   summary(bayes_multilevel)
+```{r fitrobu, class.source='klippy'}
+   mod <- robumeta::robu(formula=yi~1, data=data, studynum=study_id, var.eff.size=vi, method = "HIER", small = FALSE)
+   print(mod)
 ```
 
-As expected, our credible intervals get a little bit wider. 
+With this simple (and well balanced) data, our RVE approaches don't change the results much, but this won't always be the case. 
 
-# Correction by Modeling the Entire Sampling Covariance Matrix
+# Correction 2:  Modeling the Entire Sampling Covariance Matrix
 
-Of course, we can also take an approach proposed by @Noble2017, where we fit the covariance matrix directly by simply assuming that effects that come from the same study are correletd by r = 0.5. Ultimately, one could change this correlation, depending on the situation and context, but r = 0.5 will probably suffice in many situations. This assumes, however, that the degree of correlation among effect sizes within a study is the same across studies. This assumption is relaxed in the RVE approaches decsribed above. We can also test whether this is a safe assumption by combining it with a ClubSandwich estimator. We can build the matrix an implement this approach as follows:
+Of course, we can also take an approach proposed by @Noble2017, where we fit the covariance matrix directly by simply assuming that effects that come from the same study are correlated by r = 0.5. Ultimately, one could change this correlation, depending on the situation and context, but r = 0.5 will probably suffice in many situations. This assumes, however, that the degree of correlation among effect sizes within a study is the same across studies. This assumption is relaxed in the RVE approaches described above. We can also test whether this is a safe assumption by combining it with a ClubSandwich estimator. We can build the matrix an implement this approach as follows:
 
-```{r VCVmatrix, class.source='klippy'}
+```{r VCVmatrix, class.source="klippy"}
                  vcv <- impute_covariance_matrix(vi = data$vi, cluster = data$study_id, r = 0.5)
   
   mod_multilevel_vcv <- metafor::rma.mv(yi=yi, V = vcv, mods=~1, random=list(~1|study_id,~1|obs), data=data, test="t")
@@ -168,6 +153,20 @@ Of course, we can also take an approach proposed by @Noble2017, where we fit the
   print(mod_multilevel_vcv)
 ```
 
+# Correction 3: Applying Bayesian Multi-level Meta-analytic Model
+
+As we describe in our comment, Bayesian approaches, assuming one has a good sample size, do a very good job correcting for inflated Type I errors across a variety of situations. Bayesian MLMA models can be fit in various packages. Probably the most flexible for meta-analyst's are `MCMCglmm` [@Hadfield2010] and `brms` [@Brkner2017; @Brkner2018]. Here, we demonstrate how to fit the same MLMA model using `MCMCglmm` which has a syntax that is different from the typical one used in packages such as `metafor` and `lme4` [@Bates2015], which meta-analysts might be more accustomed too. 
+
+```{r bayes, class.source="klippy"}
+    prior <- list(R = list(V = 1, nu = 0.002), 
+                  G = list(G1 = list(V = 1 , nu = 1, alpha.mu=0, alpha.V=25^2)))
+
+   bayes_multilevel <- MCMCglmm(yi ~ 1, mev = data$vi, random = ~ study_id, data = data, prior = prior, verbose = FALSE)
+   summary(bayes_multilevel)
+```
+
+As expected, our credible intervals get a little bit wider. Bayesian models are the most conservative here given this is not a large data set.
+
 # Conclusions
 The goal of our short tutorial was to dispel the idea that overcoming, and implementing, solutions to deal with non-independent effect sizes when working with multi-level meta-analytic models is challenging. Our simulations show [@Nakagawa2021] that there are a number of very easily implemented solutions. As such, meta-analyst's can harness the power of MLMA models without the need to average effect sizes, as suggested by @Song2020.