---
title: "Dropout Rate Meta-analysis of Substance Abuse Studies"
author:
- Sara Lappan, PhD^[University of Alabama at Birmingham]
- '**Code:** Andrew W Brown, PhD^[Indiana University School of Public Health-Bloomington]'
- Peter Hendricks, PhD^[University of Alabama at Birmingham]
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_document:
    number_sections: yes
    toc: yes
    toc_depth: 4
    toc_float:
      collapsed: no
  word_document:
    toc: yes
    toc_depth: '4'
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE
                      , eval=TRUE
                      , comment = ""
                      , include=TRUE
                      , fig.height = 7)

long <- TRUE #change output for whether to have long or short report
##Set parameters
profile <- long  #whether or not to do profile plots
plots <- TRUE  #whether or not to plot meta-regressions and sub-group meta-analyses

```

# Information
This is the full data analysis code for the Dropout Meta-analysis.

**Note**: *Confirming individual pieces of this code is challenging if running interactively. If the `dropout.Rmd` file is opened and `knit` in its entirety there have been no problems. When the code is run chunk by chunk, on the other hand, the call to `knit_child` throws an error because it is apparently not designed to be run interactively. The environment successfully stores the information needed for the child document to run, so this error can be circumvented by expanding the `for` loops for each variable and running the child code manually if need be.*

```{r set environment, include=long, results='asis'}
library(xlsx) #to be able to import Excel file
library(metafor) #package for meta-analyses
library(knitr) #combining code and markdown into HTML
library(pander) #adding text formatting to chunks
library(kableExtra) #additional table formatting
source("functions.R") #source the custom functions for this project
source("datadictionary.R") #interpretation of variables to make more readable output

```

```{r echo=FALSE, include=long, results='asis'}
pandoc.header("Session Info", level=2)
```

```{r sessioInfo, echo=FALSE, include=long}
print(sessionInfo())
```

# Data import and cleaning
Data were imported from Excel, extraneous columns were removed, and incidental empty rows and columns removed.

```{r import and format, echo=long, include=long}
#Import data from Excel
filename <- "Meta Dataset 190530.xlsx" #the filename is manually called by date

#import dropout data from file
main.df <- read.xlsx(file = filename,
                     sheetName = "Study-level data")

#import the substance use data from file
subx.df <- read.xlsx(file = filename,
                     sheetName = "Frequency and Quantity")
#import the substance use data from file
diag.df <- read.xlsx(file = filename,
                     sheetName = "Comorbid Diagnoses")

#clean off extraneous columns
extra.columns <- c("Income_.annual","ASI","Country")

main.df <- main.df[,-which(colnames(main.df) %in% extra.columns)]

#clean off NA rows; this assumes there is always an identifier in the first column
dataRow <- !is.na(main.df[,1])
main.df <- main.df[dataRow,]

dataRow <- !is.na(subx.df[,1])
subx.df <- subx.df[dataRow,]

dataRow <- !is.na(diag.df[,1])
diag.df <- diag.df[dataRow,]

#Fix the name of 'Arm #'
colnames(main.df)[which(colnames(main.df) == "Arm..")] <- "Arm" #rename
```

## Cleaning and checking study-level variables
After data were imported and the data frame formatted, several important checks were conducted.

First, each study should have a unique ID within a paper. If non-unique arms were detected, they were flagged for evaluation.

```{r study-level cleaning, echo=long, include=long}
#Check for multiple of the same arm#
repArms <- NULL
for(i in unique(main.df$ID)){
  tmp <- main.df$Arm[which(main.df$ID == i)]
  if(length(unique(tmp)) != length(tmp)) #each arm ID should be unique within paper ID
    repArms <- c(repArms, as.character(main.df$Author[which(main.df$ID == i)])[1])
}
if(length(repArms) > 0) warning(paste("Stud(y/ies) with repeated study arm designations:"), repArms)
#The lack of a warning indicates that all is well.  
```

Next, the proportion of participants dropping out is checked. The data were extracted as a proportion, so any value less than 0 or greater than 1 is impossible.

```{r dropout proportion, echo=long, include=long}
#any greater than 1 have to be percentages, which indicates an error in the Excel sheet
if(any(main.df$Dropout_prop > 1)) 
  warning(paste("Some have dropout proportion > 1. Impossible. Could be entered as %. See rows",
                paste(which(main.df$Dropout_prop > 1))))
```

Because meta-analysis functions are dependent on `n` being an integer, the number who dropped out was calculated, but also checked for whether there were more estimated to have dropped out than were in the study arm total, which is impossible.

```{r dropout check, echo=long, include=long}
#calculate n dropped
nDrop <- main.df$N*main.df$Dropout_prop

#check that none somehow end up with more dropped than total
if(any(nDrop > main.df$N)) warning("Some studies have more dropped participants than total.") #should be FALSE

#check that any are less than 1 person dropped
if(any(nDrop < 1)) warning(paste("Some arms have less than 1 dropout.\n",
                                 "See papers:", paste(unique(main.df$Author[which(nDrop < 1)]), collapse = "\n")))
```

Papers returned by a warning indicated that some studies showed that less than 1 person dropped out (note that a fraction of a person could be estimated to have dropped out at this point because the number of dropouts was rounded to integers later). Papers with such low dropouts were checked and confirmed.

```{r round dropouts, echo=long, include=long}
main.df$Dropout_n_raw <- nDrop #do not round nDrop
main.df$Dropout_n <- round(nDrop)
```

##Data inclusion/exclusion rules

**Nine data inclusion/exclusion cases were identified.**

|Case|Dropout|Moderator|Decision/Notes|
|--|----------------|-----------------|----------------------------------------|
|1|Averaged over included and excluded arms|All cases|Exclude Study<br> Dropout 'tainted' by excluded groups|
|2|Reported separately per arm|None|Include as arm nested in study in hierarchical analysis|
|3|Pooled only over included arms|None|Include as single value for study (remove replicate arms)|
|4|Reported separately per arm|Reported separately per arm|Include as arm nested in study in hierarchical analysis|
|5|Pooled only over included arms|Reported separately per arm; non-poolable (categorical) moderator|Exclude Study<br> Cannot pool across separate categories (e.g., average over male/female)<br>  Including it in both categories separately would bias toward the null and violate the independence assumption|
|6|Pooled only over included arms|Reported separately per arm; Poolable (continuous) moderator|Include after pooling moderator<br> Assumes the average 'effect' of two moderator values would converge on the pooled dropout(e.g., a single pooled dropout rate for multiple moderator values)<br> Unclear if this induces regression dilution/attenuation bias|Similar to case 8|
|7|Reported separately per arm|Pooled across arms|Include as arm nested in study in hierarchical analysis <br> Analogous to school being treated with an intervention and measuring the outcome on multiple classrooms. Not all studies have this structure, and in some cases very few are pooled<br> Unclear if this inflates variance (e.g., multiple dropout rates for a single pooled moderator value)|
|8|Pooled only over included arms|Pooled only over included arms|Include single value for study|
|9|Any|Multiple values for the same comorbid diagnosis or substance class|Exclude Arm<br> Example: a study reports proportion of participants using two different opioids or with two mental comorbid diagnoses separately. There is no way to pool the values because the proportions may refer to the same or different individuals.|

The following processing accounts for these various cases. The `plo.rma.mv` function also accounts for some of these pooling rules dependent on whether the moderator is continuous or categorical through the `pool.cleaning` function.

##Processing averaged dropout rates
Some studies reported dropout rates averaged across multiple study arms. In some analyses, this can still be used. However, if a study involved study arms that were excluded from our analysis and pooled the dropout rate, then the study has to be excluded. There is no way to get the dropout rate only for arms we are interested in. This is most easily detected if an average dropout rate is reported but only one study arm is included. (See Case 1 above).

```{r pooled dropout check, echo=long, include=long}
#get study IDs for studies with total dropout reported
byStudy <- unique(main.df$ID[which(main.df$Dropout_Rep == 0)]) 

#create new instance of main.df to pool dropouts with arms only reported by study
clean.df <- main.df

#remove single arm studies with average dropouts reported
one.arm <- NULL
for(i in byStudy){
  arms <- which(clean.df$ID == i) #find arms with a given ID
  if(length(arms) == 1) {
    one.arm <- c(one.arm, as.character(clean.df$Author[which(clean.df$ID == i)]))
    clean.df <- clean.df[-which(clean.df$ID == i),]
  }
}

if(!is.null(one.arm)) 
  warning(paste0("'Dropout averages' requires all arms to be included from paper. ",
                 "The following study was REMOVED: ",
                 one.arm))

clean.df <- clean.df[-which(clean.df$Author == "Stevens & Hollis (1989)"),]

```

All studies with average dropout rates and multiple study arms were manually checked to be sure that all study arms in the original paper were eligible for inclusion for this analysis. Only Stevens & Hollis (1989) was removed for this reason.

```{r manual checking, echo=long, include=long}
#use the following code to get other study names for manual checking
#get all studies reported pooled dropout
tmp <- unique(main.df$Author[which(main.df$ID %in% byStudy)])

#only get studies that are not one-arm studies, which are excluded, or the "Stevens & Hollis", which is established as excluded
tmp <- tmp[-which(tmp %in% c(one.arm, "Stevens & Hollis (1989)"))]
```


## Cleaning and formatting study-arm-level substance and comorbid diagnosis data
Data from the substances and comorbid diagnoses worksheets were merged with the study-level data by first confirming that the Author (year) lists matched perfectly.

```{r study-arm cleaning, echo=long, include=long}
#Rename first two columns for consistency
colnames(subx.df)[c(1,2)] <- c("Author", "Arm")
colnames(diag.df)[c(1,2)] <- c("Author", "Arm")

#remove studies not in the 'clean.df' set
subx.df <- subx.df[which(subx.df$Author %in% unique(clean.df$Author)),]
diag.df <- diag.df[which(diag.df$Author %in% unique(clean.df$Author)),]

###merge###
#check that the author lists perfectly overlap
if(length(setdiff(subx.df$Author, clean.df$Author)) > 0) warning("Author (year) lists do not match between clean.df and subx.df.")
if(length(setdiff(diag.df$Author, clean.df$Author)) > 0) warning("Author (year) lists do not match between clean.df and diag.df.")

subclean.df <- merge(x = clean.df, 
                     y = subx.df, 
                     by = c("Author","Arm"),
                     all = T)

diagclean.df <- merge(x = clean.df, 
                      y = diag.df, 
                      by = c("Author","Arm"),
                      all = T)

```

Each study arm in the merged set should have at least one substance, given that all studies were about substance abuse. On the other hand, not all studies will have reported comorbid diagnoses, so those arms were removed from the merged comorbid diagnosis data set. Typically, if one study arm reported a comorbid diagnosis the others should, too; if they did not, they were flagged for follow-up.

```{r study-arm cleaning2, echo=long, include=long}
#Check whether all arms have a substance listed (they should)
emptysub <- which(is.na(subclean.df$Substance))
if(length(emptysub) > 0){
  warning(paste("The following have no substances listed:", 
                paste(unique(subclean.df$Author[emptysub]), collapse = ", ")))
  
  #remove empty substances to keep the analysis moving
  subclean.df <- subclean.df[-emptysub,]
}

#Find arms with no comorbid diagnoses; it is expected some will not
emptydiag <- which(is.na(diagclean.df$DX))
emptydiag.authors <- unique(diagclean.df$Author[emptydiag]) #get authors

#check for articles where all arms do not have diagnoses; lack of diagnoses could be a problem
diag.mix <- NULL #object to hold author list for studies with only some diagnoses
for(a in emptydiag.authors){
  tmp <- diagclean.df$DX[which(diagclean.df$Author == a)] #get diagnoses for each article
  tmp2 <- is.na(tmp) #determine if they are NA
  if(sum(tmp2) < length(tmp)) diag.mix <- c(diag.mix, a) #if there are fewer NAs than total rows, then it's mixed
}

if(length(diag.mix) > 0){
  warning(paste("The following have arms with and without comorbid diagnoses:", 
                paste(diag.mix, collapse = ", ")))
}

if(length(emptydiag) > 0){
  warning(paste("The following studies have from one to all arms without a comorbid diagnosis:\n", 
                paste(emptydiag.authors, collapse = "\n")))
  
  #remove empty diagnoses; this removes the diag.mix, too
  diagclean.df <- diagclean.df[-emptydiag,]
}
```

Occasionally, studies will report more than one substance in the same category within a study arm, such as multiple opiates. Unfortunately, the values are not able to be included easily into models because of dependence of observations (i.e., the same subjects would be counted twice for the same substance; see Case 9 above). The overlap can also be because of a data entry error, so they were flagged and confirmed.

Similarly, multiple of the same category of comorbid diagnoses may be listed, so they were also confirmed.

```{r study-arm cleaning3, echo=long, include=long}
##Check for the same substance being listed more than once within arm##
multisub <- NULL #object to hold info on multiple instances of a substance within arm
for(auth in unique(subclean.df$Author)){ #go through each study
  tmp <- subclean.df[which(subclean.df$Author == auth),] #get all arms
  test <- (length(unique(tmp$Arm)) * length(unique(tmp$Substance))) == length(tmp$Arm) #test for multisub
  if(!test) multisub <- c(multisub, as.character(tmp$Author)[1]) #if so, add study to list
}
#View(subclean.df[which(subclean.df$Author %in% multisub),c(1,2,3,70,71,72,73,74,75,76)])
#NOTE: the following are confirmed to have multiple of same substance category, like different opiates
# "Coviello et al. (2006)"
# "Iguchi et al. (1997)"
# "Morgenstern et al. (2009)"

# Show which diagnoses in the above three papers are replicated within arms
# sapply(X = multisub, FUN = function(X) table(subclean.df$Substance[which(subclean.df$Author == X)]))
# #Coviello: substance 4
# #Iguchi: substance 4
# #Morgenstern: substance 10

##Check for the same diagnosis being listed more than once within arm##
multidiag <- NULL #object to hold info on multiple instances of a diagnosis within arm
for(auth in unique(diagclean.df$Author)){ #for each study
  tmp <- diagclean.df[which(diagclean.df$Author == auth),] #get diagnoses
  test <- (length(unique(tmp$Arm)) * length(unique(tmp$DX))) == length(tmp$Arm) #test for multidiag
  if(!test) multidiag <- c(multidiag, as.character(tmp$Author)[1]) #if so, add study to list
}
#View(diagclean.df[which(diagclean.df$Author %in% multidiag),c(1,2,3,70,71,72,74)])
#Confirmed multiple diagnoses to be expected

##Ordering data sets
clean.df <- clean.df[order(clean.df$Dropout_prop),] # sort by dropout rate
subclean.df <- subclean.df[order(subclean.df$Dropout_prop),] #sort by dropout rate
diagclean.df <- diagclean.df[order(diagclean.df$Dropout_prop),] #sort by dropout rate

```

##Final counts of studies included
The PRISMA diagram has information regarding inclusion/exclusion from the point of search. Below is inclusion/exclusion summaries starting from the final data file.

* **Studies in imported file**: `r length(unique(main.df$ID))` 
* **Study arms in imported file**: `r nrow(main.df)` 
* **Studies in processed data set**: `r length(unique(clean.df$ID))` 
* **Study arms in processed data set**: `r nrow(clean.df)` 
* **Dropout rates in processed data set**: `r nrow(pool.cleaning(clean.df, "ID"))` 
* **Total subjects in processed data set**: `r fmt0(sum(pool.cleaning(clean.df, "ID")$N))` 

**NOTE:** 

* The difference between study arms and dropout rates is because some studies include only pooled dropout rates over multiple study arms.
* The number of included studies and dropout rates varies from analysis to analysis dependent on whether the moderator was reported in a given study. The number of included studies and dropout rates are reported for each analysis.

#Estimating overall dropout rates
Several estimates were made of the overall dropout rate, from the simplistic to the appropriately complex. These multiple models were fit to demonstrate the variability from one method to the next, but also to demonstrate that the dropout rates among methods are not unreasonably different from each other.

##Manual fixed-effects meta-analysis of proportions
The overall dropout rate was manually calculated as a fixed-effects, individual-value meta-analysis, not taking into consideration that study arms are related within study for studies that reported dropout rates for each arm independently. This was accomplished by simply looking at the proportion of dropouts summed across all studies divided by the total number of participants across all studies. This method makes three unreasonable assumptions: there is no random effect of study; that study arms are not clustered within study; and (taken together) that each individual across studies essentially is a completely random and independent sample from the population.

```{r manual proportion, echo=long}
##OVERALL MANUAL CALCULATION
df <- pool.cleaning(clean.df, "ID") #adjust all avg dropouts
n <- sum(df$N) #get total number of participants
ns <- sum(df$Dropout_n_raw) #get total number that dropped out
nf <- n - ns #get number who did not dropout
z <- qnorm(0.975) #set z to threshold
pred <- ns/n #calculate proportion dropped out
ci <- z*sqrt((1/n)*pred*(1-pred)) #calculate the confidence interval
ci.lb <- pred - ci #calculate the lower confidence interval
ci.ub <- pred + ci #calculate the higher confidence interval
manual.raw <- c(pred, ci.lb, ci.ub) #create the output object
names(manual.raw) <- c("Dropout","Lower 95% CI","Upper 95% CI") #name the output object
manual.round <- round(manual.raw, 4) #round the object for display
manual.round <- manual.round * 100 #convert to percent
cat("\n")
kable(as.data.frame(t(manual.round)))%>%
  kable_styling(full_width = F)
```

##Manual, inverse variance, fixed-effects, logit meta-analysis 
The overall dropout rate was then calculated by first estimating the logit proportion and variance, followed by manual calculation of the mean effect size (proportion dropout) and its variance. This approach makes the untested assumption that there is no heterogeneity among studies, and that study arms are independent (i.e., ignores any ICC among study arms within studies).

```{r manual inverse variance fixed, echo=long}
#calculate the logit proportion dropout and its variance for each dropout rate
es <- escalc(measure = "PLO",
             xi = df$Dropout_n,
             ni = df$N)

es <- as.data.frame(es) #convert to data.frame for ease of use

dat <- cbind(rownames(df),as.data.frame(es), df) #add the rownames as a column and append the original df
colnames(dat)[1] <- "arm"

mean <- sum(es[1]/es[2])/sum(1/es[2]) #calculate the overall mean effect size
ci <- sqrt(1/sum(1/es[2]))*z #calculate the 95% interval on the logit scale
iv.raw <- c(mean, mean-ci, mean+ci) #create an object with the mean and 95% CI on logit scale
iv.raw <- transf.ilogit(iv.raw) #transform to proportion
names(iv.raw) <- c("Dropout", "Lower 95% CI", "Upper 95% CI") #name the columns
iv.round <- round(iv.raw, 4) #round the object for display
iv.round <- iv.round * 100 #convert to percent

kable(t(iv.round))%>%
  kable_styling(full_width = F)
```

##Overall dropout by fixed-effects meta-analysis using `rma.mv`
Next, the overall dropout rate was calculated using the `rma.mv` function of the `metafor` package. The same assumptions are made as in the above manual calculation.

```{r overall fe rma.mv, echo=long}
###OVERALL RMA.MV###
#use df and dat from above to calculate fixed-effects, inverse variance meta-analysis
out <- rma.mv(yi, 
              vi, 
              data = dat)
print(out)
rma.fe.out <- predict(out, transf=transf.ilogit) #convert to proportion
rma.fe.out <- as.numeric(unlist(rma.fe.out)[c(1,3,4)]) * 100 #convert to percent
rma.fe.out <- round(rma.fe.out,2)
names(rma.fe.out) <- c("Dropout", "Lower 95% CI", "Upper 95% CI")
kable(t(rma.fe.out))%>%
  kable_styling(full_width = F)
```

##Overall random-effects
Building onto the fixed-effects models, a random-effects model was fit using `rma.mv` that continues to assume each dropout rate is independent. This also applies to studies in which the dropout rate was pooled.

```{r overall random-effects, echo=long}
###OVERALL random-effects RMA.MV###
#use df and dat from above to calculate random-effects, inverse variance meta-analysis
out <- rma.mv(yi, 
              vi, 
              random = ~1 | arm,
              data = dat)
print(out)
rma.re.out <- predict(out, transf=transf.ilogit) #convert to proportion
rma.re.out <- as.numeric(unlist(rma.re.out)[c(1,3,4,5,6)]) * 100 #convert to percent
rma.re.out <- round(rma.re.out,2)
names(rma.re.out) <- c("Dropout", "Lower 95% CI", "Upper 95% CI", "Lower 95% PI", "Upper 95% PI")
kable(t(rma.re.out))%>%
  kable_styling(full_width = F)
```

##Overall dropout by multi-level, random-effects meta-analysis
Finally, the overall dropout rate was estimated using a multilevel random-effects meta-analysis using a logit transformed proportion. This was accomplished through the custom `plo.rma.mv` wrapper function for the `rma.mv` function of `metafor` (see `functions.r` file). This model also included the hierarchical nature of the data, with study arms being nested with studies. Note that in cases where dropout rates were pooled, the study was reported as a single entry and thus no within-study correlation could be calculated.

```{r overall rma.mv, echo=long}
overall.re.out <- plo.rma.mv(clean.df
                             ,varname = NULL
                             ,title = "Overall",
                             profile = profile,
                             long = long)
print(overall.re.out$model)
cat(paste("\nICC:", fmt4(overall.re.out$ICC)))
cat(paste("\nI^2.1", fmt1(overall.re.out$I2[1])))
cat(paste("\nI^2.2", fmt1(overall.re.out$I2[2])))
kable(overall.re.out$prediction)%>%
  kable_styling(full_width = F)
if(plots){
  forest(
    x = overall.re.out$model$yi[order(overall.re.out$model$yi)] #sort by point estimate
    ,vi = overall.re.out$model$vi[order(overall.re.out$model$yi)] #sort by point estimate
    ,annotate = FALSE
    ,at=c(0,0.2,0.4,0.6,0.8,1)
    ,transf = transf.ilogit
    ,xlim = c(-0.1, 1.1)
    ,alim = c(0,1)
    ,refline = overall.re.out$prediction$Dropout/100
    ,slab = NA
    ,xlab = "Dropout (proportion)"
    ,cex.lab = 1
    ,cex.axis = 1
    ,psize = 10
  )
  
  #add prediction interval under diamond
  lines(c(overall.re.out$prediction$`Lower 95% PI`/100,
          overall.re.out$prediction$`Upper 95% PI`/100)
        ,c(-5,-5))  
  
  #manually add the summary diamond
  addpoly(x = overall.re.out$prediction$Dropout/100
          ,ci.lb = overall.re.out$prediction$`Lower 95% CI`/100
          ,ci.ub = overall.re.out$prediction$`Upper 95% CI`/100
          ,rows = -5
          ,efac = 1
          ,cex=1.2
          ,col='red'
          ,annotate=FALSE
  )
  
  #Add annotation
  text(x=c(0,0,0)
       ,y=c(25,15,5)
       ,labels = c(
         paste0("Mean Dropout: "
                ,fmt2(overall.re.out$prediction$Dropout/100)
         ),
         paste0("95% CI: "
                ,fmt2(overall.re.out$prediction$`Lower 95% CI`/100)
                ,","
                ,fmt2(overall.re.out$prediction$`Upper 95% CI`/100)
         ),
         paste0("95% PI: "
                ,fmt2(overall.re.out$prediction$`Lower 95% PI`/100)
                ,","
                ,fmt2(overall.re.out$prediction$`Upper 95% PI`/100)
         )
       )
       ,adj=0
  )
  
}
```

```{r overall meta-regression legend, echo=FALSE, results='asis'}
if(plots) cat("\n**Figure:** Overall, random-effects, hierarchical, logit meta-analysis. Each horizontal line represents the point estimate and 95% CI of each study, back-transformed from logit scale. Symbols demark the point estimate, and are not proportional to study or study-arm weight. The dashed line represents the mean proportion dropout across studies. The polygon represents the average proportion dropout and its 95% CI, with a 95% prediction interval.")
```

##Conclusion about overall dropout rate
The manual and random-effects analyses are not identical, as would be expected. Although there are differences among them, they are not large differences, giving some degree of face validity. The model with the highest theoretical support is the hierarchical random-effects model, and is considered the primary analysis for the overall dropout rate.

Importantly, there remains significant residual heterogeneity, as is expected given the wide variety of studies included. To investigate sources of this heterogeneity, a series of moderator meta-regressions (for continuous predictors) or sub-group meta-analyses (for categorical preditors) are described below.

The interpretation of the dropout rates and 95% CI are limited in the random-effects models, as indicated by the 95% Prediction Interval (PI). The sizeable heterogeneity indicates that although the average across all studies is expected to be around `r round(overall.re.out$prediction[1],0)`%, the expectation for any given study is much less predictable, with a 95% range of study dropout rates from about `r round(overall.re.out$prediction[4],0)`% to `r round(overall.re.out$prediction[5],0)`%.

#Testing study-level moderators
To evaluate whether characteristics of the studies could help explain the heterogeneity scene in the overall models, moderator sub-group meta-analyses and meta-regressions were calculated. This section describes study-level moderators, which are moderator variables that were constant at the study level, meaning they did not vary from study-arm to study-arm.

##Categorical study-level moderators
Sub-group meta-analyses were conducted for moderator variables that were categorical. These models used the same inverse-variance, hierarchical approach as the overall model, but included categorical, study-level variables for sub-group analyses. For each moderator analysis, profile plots may be created to evaluate the two variance components (within and between study), and the intraclass correlation (ICC) is also reported. In cases where the ICC is 0 or 1, the profile plots for one of the variance components is maxed out at 0.

```{r categorical-moderators, echo=long}
###Test moderators###
#Moderators from worksheet 1 only
###categorical variables###
catvars <- c(
  "Inter_degree",
  "Type",
  "Pregnant",
  "Manualized",
  "Setting",
  "Pharmacotherapy",
  "TX_orientation",
  "Limited",
  "Training.for.fidelity",
  "Format",
  "Eff",
  "Dependence",
  "Treatment.Seeking",
  "Country_Dev"
)
out <- NULL
categorical.table <- NULL #instantiate table to be printed at end
for(v in catvars){
  outname <- common_varname[match(v,common_varname),2] #get readable version of varname
  if(is.na(outname)) outname <- v #if it does not exist, revert to variable name
  
  cat.out <- suppressMessages(plo.rma.mv(df = clean.df,
                                         varname = v,
                                         continuous = FALSE,
                                         title = outname,
                                         profile = profile,
                                         long = long)
  )
  assign(paste0(v,".out"), cat.out)
  out <- c(out, knit_child('categorical.Rmd', quiet=TRUE))
  categorical.table <- rbind(categorical.table, cat.table.out)
}

```

```{r print categorical, echo=FALSE, results='asis'}
cat(paste(out, collapse = '\n'))
```

## Continuous study-level moderators
Continuous study-level moderators were evaluated using meta-regressions, again based on the inverse-variance, hierarchical, random-effects approach described above. Profile plots were also fit as described in the 'Categorical' section, and meta-regression plots were created.

```{r continuous-moderators, echo=long}
##continuous variables##
#can add predvals, but presently just uses the six categories from summary() function
contvars <- c(
  "Age",
  "Gender_Male",
  "Education_years",
  "White",
  "Hispanic",
  "AA",
  "Race_eth_other",
  "Adjusted_income_annual",
  "Marital_not_married",
  "Unemployment",
  "Employment_other",
  "Inter_exp_years",
  "Year",
  "Sessions",
  "Sess_length_min",
  "Session_period_weeks"
)
out<-NULL
continuous.table <- NULL #instantiate table to be printed at end
for(v in contvars){
  outname <- common_varname[match(v,common_varname),2] #get readable version of varname
  if(is.na(outname)) outname <- v #if it does not exist, revert to variable name
  
  cont.out <- suppressMessages(plo.rma.mv(df = clean.df,
                                          varname = v,
                                          continuous = TRUE,
                                          title = outname,
                                          profile = profile,
                                          long = long)
  )
  assign(paste0(v,".out"), cont.out)
  out <- c(out, knit_child('continuous.Rmd', quiet=TRUE))
  continuous.table <- rbind(continuous.table, cont.table.out)
  
}
```

```{r print continuous, echo=FALSE, results='asis'}
cat(paste(out, collapse = '\n'))

```

#Study-arm level moderators: Substances
Data were also extracted on study-arm level moderators. This included characterizing the participants' target substance use, off-target substance use, and comorbid diagnoses. In this section, we desscribe the substances. 

Simultaneous moderation of substances is not possible with these data because of independence issues, namely that the same study arm could have reported more than one substance and would therefore be counted twice. We chose to test each substance separately, excluding arms with more than one substance in the same class (e.g., multiple different opiates). Continuous and categorical moderators were fit as described above.

```{r substances, echo=long}
out<-NULL #placeholder object for knitr output
sx.table <- NULL #placeholder object for substance tables
sub.deg.table <- NULL #placeholder object for substance degree of use tables

for(sx in sort(unique(subclean.df$Substance))){ #for each substance
  substance.df <- subclean.df[which(subclean.df$Substance == sx),] #extract only ones with substance sx
  
  ##exclude arms with multiple substances per arm
  exclude <- NULL #new object to get exclusion id/arm pairs; column 1 is substance, 2 is arm
  for(s in unique(substance.df$ID)){ #for each study
    tmp <- substance.df[which(substance.df$ID == s),] #get one study at a time
    for(a in unique(tmp$Arm)){ #for each arm within study
      tmp2 <- tmp[which(tmp$Arm == a),] #get a single arm in the study
      if(nrow(tmp2) != 1) { #if that arm has more than one instance of a substance
        exclude <- rbind(exclude, c(s,a)) #add to exclusion list
      }
    }
  }
  
  if(!is.null(nrow(exclude))){ #if there is an exclusion list
    for(e in 1:nrow(exclude)){
      substance.df <- substance.df[-which(substance.df$ID == exclude[e,1] & substance.df$Arm == exclude[e,2]),] #remove all excluded study arms
      message(paste("Arm", exclude[e,2], "of ID", exclude[e,1], "removed for multiple subtances within the class of Substance", sx))
    }
    for(id in unique(exclude[,1])){ #re-loop through to check for now-invalid pooled dropout rates
      dropout.rep <- substance.df$Dropout_Rep[which(substance.df$ID == id)] == 0 #look for pooled rates
      if(length(dropout.rep)>0) #if the id is stil in the set
        if(dropout.rep == TRUE) #and if the dropout rate is pooled
          substance.df <- substance.df[-which(substance.df$ID == id),] #remove the id
        if(long) message(paste("ID", id, "removed because of removing arms for multiple substances within the class of Substance",sx))
    }
  }
  out <- c(out, knit_child('substances.Rmd', quiet=TRUE))
}
```

```{r print substances, echo=FALSE, results='asis'}
cat(paste(out, collapse = '\n'))

```

#Study-arm level moderators: Comorbid Diagnoses
Similar to the substances, simultaneous moderation of comorbid diagnoses is also not possible with these data because of independence issues, namely that the same study arm could have more than one comorbid diagnosis and would be counted twice. Therefore, we chose to test each comorbid diagnosis separately, excluding arms with more than one comorbid diagnosis in the same class (e.g., multiple mental illnesses). Comorbid diagnoses were fit as categorical data, as described above.

```{r diagnoses, echo=long}
###Testing one diagnosis at a time###

out<-NULL #placeholder object for knitr output
diag.table <- NULL #placeholder object for diagnosis tables

for(d in 1:6){
  diagnosis.df <- diagclean.df[which(diagclean.df$DX == d),] #extract only ones with diagnosis d
  
  exclude <- NULL #new object to get exclusion id/arm pairs
  
  for(s in unique(diagnosis.df$ID)){ #for each study
    tmp <- diagnosis.df[which(diagnosis.df$ID == s),] #get one study at a time
    for(a in unique(tmp$Arm)){ #for each arm within study
      tmp2 <- tmp[which(tmp$Arm == a),] #get a single arm in the study
      if(nrow(tmp2) != 1) { #if that arm has more than one instance of a diagnosis
        exclude <- rbind(exclude, c(s,a)) #add to exclusion list
      }
    }
  }
  
  if(!is.null(nrow(exclude))){
    for(e in 1:nrow(exclude)){
      diagnosis.df <- diagnosis.df[-which(diagnosis.df$ID == exclude[e,1] & diagnosis.df$Arm == exclude[e,2]),]
    }
  }
  
  outname <- DX_dic[match(d,DX_dic),2] #get readable version of varname
  if(is.na(outname)) outname <- d #if it does not exist, revert to variable name
  #diagnosis <- paste0("**Comorbid Diagnosis: ", outname, "**") #bolded version
  diagnosis <- outname
  
  out <- c(out, knit_child('diagnoses.Rmd', quiet=TRUE))
}
```

```{r print diagnoses, echo=FALSE, results='asis'}
cat(paste(out, collapse = '\n'))
```

#Summary tables for publication
The following tables summarize the above output, divided into groups as analyzed above.

##Study-level moderator meta-regression results

```{r study-level continuous combined table, echo=long, results="asis"}
colnames(continuous.table) <- c('Moderator[note]'
                                ,'Studies (dropout rates)[note]'
                                ,'Logit regression[note]'
                                ,'Dropout at Q1[note]'
                                ,'Dropout at Median'
                                ,'Dropout at Q3'
                                ,'Residual Heterogeneity % (p)'
                                ,'Test for Moderation (p)')
kable(continuous.table) %>%
  add_header_above(header = c("Study-level moderator meta-regression results"=8), align = "l", bold=TRUE) %>%
  add_footnote(c("Separate meta-regressions were calculated for each continuous moderator.",
                 "The number of studies and dropout rates included varies by moderator because of differences in reporting across studies. When dropout rates were averaged across study-arms, either because of pooled moderators or dropout rates, the study has only a single dropout rate, regardless of how many study-arms were included.",
                 "Because the regression was calculated on the logit scale, the parameter estimates are reported on the logit scale. To use the regression to get proportion dropout, enter the value of the moderator of interest into the regression equation, and back-transform to a proportion using a function like the `transf.ilogit` function in the `metafor` package in `R`. The regression represents dropout rates as a function of the study-level or study-arm-level average moderator value, not individual-participant-level data.",
                 "Q1, Median, and Q3 refer to the first quartile, median, and third quartile (respectively) of the study-level or study-arm-level values of the moderator; note these do not necessarily represent individual-participant-level data. Because the regression is on the logit scale, these columns provide interpretable, back-transformed dropout rates in percent at Q1, Median, and Q3 moderator values. Data are presented as the **bolded** moderator value, followed by the estimated dropout percentage and 95% CI.")) %>%
  kable_styling(full_width = F)
```

##Study-level sub-group meta-analysis results

```{r study-level categorical combined table, echo=long, results='asis'}
categorical.table <- as.data.frame(categorical.table)
colnames(categorical.table) <- c('Moderator'
                                 ,'Moderator Value'
                                 ,'Studies (dropout rates)'
                                 ,'Dropout (95% CI)'
                                 ,'Residual Heterogeneity % (p)'
                                 ,'Test for Moderation (p)')
rownames(categorical.table) <- c(1:nrow(categorical.table))
kable(categorical.table) %>%
  add_header_above(header = c("Study-level sub-group meta-analysis results"=6), align = "l", bold=TRUE) %>%
  add_footnote(c("Sub-group meta-analyses were calculated for categorical moderator variables.",
                 "The number of studies and dropout rates included varies by moderator because of differences in reporting across studies. When dropout rates were averaged across study-arms, either because of pooled moderators or dropout rates, the study has only a single dropout rate, regardless of how many study-arms were included.")) %>%
  kable_styling(full_width = F)

```

##Substance use meta-regression results

```{r substance combined table, echo=long, results='asis'}
colnames(sx.table) <- c('Moderator[note]'
                        ,'Studies (dropout rates)[note]'
                        ,'Logit regression[note]'
                        ,'Dropout at Q1[note]'
                        ,'Dropout at Median'
                        ,'Dropout at Q3'
                        ,'Residual Heterogeneity % (p)'
                        ,'Test for Moderation (p)')
kable(sx.table) %>%
  add_header_above(header = c("Substance use meta-regression results"=8), align = "l", bold=TRUE) %>%
  add_footnote(c("Separate meta-regressions were calculated for each substance, as well as the use pattern (e.g., Frequency of Days Used) within substance.",
                 "The number of studies and dropout rates included varies by substance and substance use pattern because of differences in reporting across studies. When dropout rates were averaged across study-arms, either because of pooled moderators or dropout rates, the study has only a single dropout rate, regardless of how many study-arms were included.",
                 "Because the regression was calculated on the logit scale, the parameter estimates are reported on the logit scale. To use the regression to get proportion dropout, enter the value of the moderator of interest into the regression equation, and back-transform to a proportion using a function like the `transf.ilogit` function in the `metafor` package in `R`. The regression represents dropout rates as a function of the study-level or study-arm-level average moderator value, not individual-participant-level data.",
                 "Q1, Median, and Q3 refer to the first quartile, median, and third quartile (respectively) of the study-level or study-arm-level values of the moderator; note these do not necessarily represent individual-participant-level data. Because the regression is on the logit scale, these columns provide interpretable, back-transformed dropout rates in percent at Q1, Median, and Q3 moderator values. Data are presented as the **bolded** moderator value, followed by the estimated dropout percentage and 95% CI.")) %>%
  kable_styling(full_width = F)
```

##Degree of Use sub-group meta-analysis results

```{r degree of use combined table, echo=long, results='asis'}
colnames(sub.deg.table) <- c('Moderator[note]'
                             ,'Moderator Value'
                             ,'Studies (dropout rates)[note]'
                             ,'Dropout (95% CI)'
                             ,'Residual Heterogeneity % (p)'
                             ,'Test for Moderation (p)')
rownames(sub.deg.table) <- c(1:nrow(sub.deg.table))
kable(sub.deg.table) %>%
  add_header_above(header = c("Degree of Use sub-group meta-analysis results"=6), align = "l", bold=TRUE) %>%
  add_footnote(c("Sub-group meta-analyses were calculated for degree of use pattern categories within each substance type.",
                 "The number of studies and study-arms included varies by substance because of the variety of substances investigated across studies. When dropout rates were averaged across study-arms, either because of pooled moderators or dropout rates, the study has only a single dropout rate, regardless of how many study-arms were included.")) %>%
  kable_styling(full_width = F)
```

##Comorbid diagnosis sub-group meta-regression results

```{r diagnosis combined table, echo=long, results='asis'}
colnames(diag.table) <- c('Moderator[note]'
                          ,'Studies (dropout rates)[note]'
                          ,'Logit regression[note]'
                          ,'Dropout at Q1[note]'
                          ,'Dropout at Median'
                          ,'Dropout at Q3'
                          ,'Residual Heterogeneity % (p)'
                          ,'Test for Moderation (p)')
kable(diag.table) %>%
  add_header_above(header = c("Comorbid diagnosis sub-group meta-analysis results"=8), align = "l", bold=TRUE) %>%
  add_footnote(c("Separate meta-regressions were calculated for each comorbid diagnosis.",
                 "The number of studies and dropout rates included varies by diagnosis because of differences in reporting across studies. When dropout rates were averaged across study-arms, either because of pooled moderators or dropout rates, the study has only a single dropout rate, regardless of how many study-arms were included.",
                 "Because the regression was calculated on the logit scale, the parameter estimates are reported on the logit scale. To use the regression to get proportion dropout, enter the value of the moderator of interest into the regression equation, and back-transform to a proportion using a function like the `transf.ilogit` function in the `metafor` package in `R`. The regression represents dropout rates as a function of the study-level or study-arm-level average moderator value, not individual-participant-level data.",
                 "Q1, Median, and Q3 refer to the first quartile, median, and third quartile (respectively) of the study-level or study-arm-level values of the moderator; note these do not necessarily represent individual-participant-level data. Because the regression is on the logit scale, these columns provide interpretable, back-transformed dropout rates in percent at Q1, Median, and Q3 moderator values. Data are presented as the **bolded** moderator value, followed by the estimated dropout percentage and 95% CI.")) %>%
  kable_styling(full_width = F)
```

