Verbose.Rmd

---
title: "Verbose Meta-Analysis Example"
bibliography: R-Pckgs.bib
csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa.csl
output:
  html_document: default
  pdf_document: default
  word_document: default
nocite: |
  @*
always_allow_html: yes
---

```{r echo=FALSE}
Author <- "Your Name"
Rev.Title <- "Review Title"
filename.ES.Data<- "alt.csv"
filename.PRISMA.Data<- "prisma.csv"
X<- "Anxiety"
Y <- "Exercise"
Z <- "Waitlist Control"


### Specify the columns containing the Study ID and means, SDs, and Ns for intervention and control within the file

col.int.means     <- "m1"
col.int.sds       <- "sd1"
col.int.ns        <- "n1"
col.cont.means    <- "m2"
col.cont.sds      <- "sd2"
col.cont.ns       <- "n2"
col.study.id      <- "StudyID"


### Specify the effect size measure

## Options in this template are SMD (standarised mean differences) and MD (raw mean differences). 
## In most cases SMD is most appropriate.

measure <- "SMD"

### Specify the model type. 
## In most cases REML should be the default

## Options in this template are 

# method <- "FE" # Fixed effect meta-analysis
# method <- "REML" # Default random effects meta-analysis
# method <- "DL" #  DerSimonian-Laird estimator
# method <- "HE" # Hedges estimator
# method <- "HS" # Hunter-Schmidt estimator
# method <- "SJ" # Sidik-Jonkman estimator
# method <- "ML" # maximum-likelihood estimator
# method <- "REML" # restricted maximum-likelihood estimator
# method <- "EB" # empirical Bayes estimator
# method <- "PM" # Paule-Mandel estimator
# method <- "GENQ" # generalized Q-statistic estimator

method <- "REML"


```


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(metafor)
library(knitr)
library(RCurl)
library(dplyr)
library(rmarkdown)
library(DiagrammeR)
#library(DiagrammeRsvg)
#library(magrittr)
#library(svglite)
#library(rsvg)
#library(png)


output<-default_output_format(knitr::current_input()) ## This will fail if run manually. Don't worry!


## These lines of code download and run the metafor_tidiers functions that implement broom type tidy data functions for rma objects
source("metafor_tidiers.R")

## This opens the PRISMA function adapted from the prismastatement package
source("prisma.R")

# Set so that long lines in R will be wrapped:
opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE)

prisma.data <- read.csv(filename.PRISMA.Data)
```

```{r calculate_ES, include=FALSE}
dat <- read.csv(filename.ES.Data, stringsAsFactors = FALSE)
  dat_ES <-
    escalc(
    measure = measure,
    m1i = get(col.int.means),
    sd1i = get(col.int.sds),  
    n1i = get(col.int.ns),  
    m2i = get(col.cont.means),
    sd2i = get(col.cont.sds),  
    n2i = get(col.cont.ns),  
    data = dat
    )

```

```{r conduct_MA, echo=FALSE}
dat_MA<- rma(yi, vi, data=dat_ES, slab=get(col.study.id))
model<-tidy.rma(dat_MA)
het.small<-glance.rma(dat_MA) %>% select(one_of(c("k", "tau2", "se.tau2", "QE", "QEp", "I2")))
eggers<-regtest(dat_MA)

```

## Data analysis

This file documents the analyses conducted for `r Author` (`r (format(Sys.Date(), "%Y"))`) *`r Rev.Title`*. Analyses were conducted using the file `r filename.ES.Data`. 

A random-effects meta-analysis (k = `r dat_MA$k`) was conducted using the `r dat_MA$method` estimator. The effect size is the standarised mean difference (hedges g).

#Results
## Studies
This review was reported according to the PRISMA statement. The flow of studies through the review is shown in the PRISMA Flowchart below (Figure 1). From the `r prisma.data$no_dupes` studies that underwent screening, a total of `r prisma.data$quantitative` were suitable for meta-analysis.


```{r prisma_png, fig.cap="Figure 1. PRISMA Flowchart", echo=FALSE, warning=FALSE, out.width="400px", message=FALSE, eval=(output$name!="html_document")}
prisma_graph<-prisma(found = prisma.data$found,
        found_other = prisma.data$found_other,
        no_dupes = prisma.data$no_dupes, 
        screened = prisma.data$screened, 
        screen_exclusions = prisma.data$screened_exclusions, 
        full_text = prisma.data$full_text,
        full_text_exclusions = prisma.data$full_text_exclusions, 
        qualitative = prisma.data$qualitative,
        quantitative = prisma.data$quantitative,
        reasons = paste(prisma.data$reasons))

# Create a PNG of this graph
export_svg(prisma_graph) %>% charToRaw %>% rsvg %>% png::writePNG('prisma_graph.png')

knitr::include_graphics('graph2.png')
```

```{r prisma_html, fig.cap="Figure 1. PRISMA Flowchart", echo=FALSE, warning=FALSE, message=FALSE, eval=(output$name=="html_document")}

prisma(found = prisma.data$found,
        found_other = prisma.data$found_other,
        no_dupes = prisma.data$no_dupes, 
        screened = prisma.data$screened, 
        screen_exclusions = prisma.data$screened_exclusions, 
        full_text = prisma.data$full_text,
        full_text_exclusions = prisma.data$full_text_exclusions, 
        qualitative = prisma.data$qualitative,
        quantitative = prisma.data$quantitative,
        reasons = paste(prisma.data$reasons))
```

## Meta-analytic synthesis

```{r forest, echo=FALSE, warning=FALSE, fig.height=het.small$k*0.5, fig.cap="Figure 2. Forest Plot"}
forest(dat_MA)
```

A random effects meta-analysis was conducted (k=`r het.small$k`) to explore the difference in `r X` between the `r Y` group and the `r Z` group. As shown in Figure 2, the average difference in `r X` was g=`r round(model$b,2)` (p=`r round(model$pval,3)`, 95% CI [`r round(model$ci.lb, 2)`, `r round(model$ci.ub, 2)`]).`r ifelse(model$pval > 0.05, "[^1]", "")` The results of this analysis are summarised in the tables below.

A Cochran's Q test was conducted to examine whether variations in the observed effect are likely to be attributable soley to sampling error (Q~(df=`r het.small$k-1`)~=`r round(het.small$QE,2)`, p=`r ifelse(het.small$QEp < 0.001, "<.001", round(het.small$QEp,3))`). `r ifelse(het.small$QEp < 0.05, "The variation in the effect is greater than would be expected from sampling error alone. It appears that the true effect varies betweeen studies.", "There is no evidence that the true effect size varies between studies.")` The I^2^ statistics indicates the proportion of variance in the observed effect attributable to sampling error. In this instance, the I^2^ = `r round(het.small$I2,2)`%.[^2]

`r ifelse(model$pval > 0.05, "[^1]: It is important to note that a p>.05 indicates lack of evidence of an effect (i.e. uncertainty) rather than evidence of no effect unless confidence intervals are sufficently narrow to rule out a clinically meaningful effect.", "")` 

[^2]: Note, this statistic is not an absolute measure of heterogeneity (although it is often interpreted as such). We strongly advise against using rules of thumb such as "small", "medium" or "large" when interpreting I^2^ values.

`r kable(cbind(model, het.small), col.names=c("*g*", "se", "z", "*p*", "95% CI LB", "95% CI UB", "k", "$\\tau$^2^", "se", "Q", "*p*", "I^2^"), row.names=FALSE, digits = 3, caption="Effect Size")`



A funnel plot was generated to allow for visual inspection of funnel plot asymmetry that can indicate reporting biases (e.g. publication bias or outcome selection bias). `r ifelse(het.small$k < 10, "However, as per Sterne and Egger (2006), statisical tests for funnel plot asymettry are underpowered when the number of studies included in the analysis is 10 or fewer. As such, a formal test of funnel plot asymmetry was not conducted", paste0("An Egger's test was conducted to detect funnel plot asymmetry (z=", round(eggers$zval, 2), ",p=", round(eggers$pval,2)))`). 


```{r fig.cap="Figure 3. Funnel Plot", echo=FALSE}
funnel<-funnel(dat_MA, back="white")
```

## Packages used in this document
```{r include=FALSE}
citPkgs <- names(sessionInfo()$otherPkgs)
write_bib(citPkgs, file="R-Pckgs.bib")
```