library(metafor) library(meta) setwd("~/Desktop/") df <- read.csv("Analysis_DepAnx.csv") names(df)[1] <- "StudyID" # Meta-analysis with only the post scores # For depression outcomes df_es_dep <- escalc("SMD", m1i = Mean_Dep_Active_End, m2i = Mean_Dep_Sham_End, sd1i = SD_Dep_Active_End, sd2i = SD_Dep_Sham_End, n1i = N_Active, n2i = N_Sham, vtype="UB", data=df) rma_dep <- rma(yi=yi, vi = vi,method= "DL", slab= StudyID ,data=df_es_dep) rma_dep forest(rma_dep, main = "Random effects meta-analysis for depression outcomes") funnel(rma_dep, main="Funnel plot for depression outcomes") regtest(rma_dep) # For anxiety outcomes df_es_Anx <- escalc("SMD", m1i = Mean_Anx_Active_End, m2i = Mean_Anx_Sham_End, sd1i = SD_Anx_Active_End, sd2i = SD_Anx_Sham_End, n1i = N_Active, n2i = N_Sham, vtype="UB", data=df) rma_Anx <- rma(yi=yi, vi = vi,method= "DL", slab= StudyID ,data=df_es_Anx) rma_Anx forest(rma_Anx, main = "Random effects meta-analysis for Anxiety outcomes") funnel(rma_Anx, main="Funnel plot for Anxiety outcomes") regtest(rma_Anx) # Doing the meta-analysis again with "meta" package rma_Dep2 <- metacont(n.e=N_Active, mean.e = Mean_Dep_Active_End ,sd.e = SD_Dep_Active_End , n.c = N_Sham, mean.c = Mean_Dep_Sham_End, sd.c = SD_Dep_Sham_End, studlab = StudyID, byvar = Target, data = df[!is.na(df$Mean_Dep_Active_End), ], sm = "SMD", method.smd= "Hedges") rma_Dep2 meta::forest(rma_Dep2, comb.fixed=FALSE) rma_Anx2 <- metacont(n.e=N_Active, mean.e = Mean_Anx_Active_End ,sd.e = SD_Anx_Active_End , n.c = N_Sham, mean.c = Mean_Anx_Sham_End, sd.c = SD_Anx_Sham_End, studlab = StudyID, byvar = Target, data = df[!is.na(df$Mean_Anx_Active_End), ], sm = "SMD", method.smd= "Hedges") rma_Anx2 meta::forest(rma_Anx2, comb.fixed=FALSE) rma_YBOCS2 <- metacont(n.e=N_Active, mean.e = Mean_YBOCS_Active_End ,sd.e = SD_YBOCS_Active_End , n.c = N_Sham, mean.c = Mean_YBOCS_Sham_End, sd.c = SD_YBOCS_Sham_End, studlab = StudyID, byvar = Target, data = df[!is.na(df$Mean_YBOCS_Active_End), ], sm = "SMD", method.smd= "Hedges") rma_YBOCS2 meta::forest(rma_YBOCS2, comb.fixed=FALSE) # Doing the meta-analysis with pre-post differences library(dplyr) df.dep <- df %>% select(StudyID, Target, Mean_Dep_Active_Baseline, SD_Dep_Active_Baseline, Mean_Dep_Active_End, SD_Dep_Active_End, N_Active, Mean_Dep_Sham_Baseline, SD_Dep_Sham_Baseline, Mean_Dep_Sham_End, SD_Dep_Sham_End, N_Sham) %>% na.exclude df.anx <- df %>% select(StudyID, Target, Mean_Anx_Active_Baseline, SD_Anx_Active_Baseline, Mean_Anx_Active_End, SD_Anx_Active_End, N_Active, Mean_Anx_Sham_Baseline, SD_Anx_Sham_Baseline, Mean_Anx_Sham_End, SD_Anx_Sham_End, N_Sham) %>% na.exclude df.YBOCS <- df %>% select(StudyID, Target, Mean_YBOCS_Active_Baseline, SD_YBOCS_Active_Baseline, Mean_YBOCS_Active_End, SD_YBOCS_Active_End, N_Active, Mean_YBOCS_Sham_Baseline, SD_YBOCS_Sham_Baseline, Mean_YBOCS_Sham_End, SD_YBOCS_Sham_End, N_Sham) %>% na.exclude c <- 0.5 # Keeping the correlation between the pre-post scores to be this, based on other studies # Calculating depression mean changes df.dep$m.dep.ch.act <- df.dep$Mean_Dep_Active_Baseline - df.dep$Mean_Dep_Active_End df.dep$sd.dep.ch.act <- sqrt(df.dep$SD_Dep_Active_Baseline^2 + df.dep$SD_Dep_Active_End^2 + 2*c*df.dep$SD_Dep_Active_Baseline*df.dep$SD_Dep_Active_End) df.dep$m.dep.ch.shm <- df.dep$Mean_Dep_Sham_Baseline - df.dep$Mean_Dep_Sham_End df.dep$sd.dep.ch.shm <-sqrt(df.dep$SD_Dep_Sham_Baseline^2 + df.dep$SD_Dep_Sham_End^2 + 2*c*df.dep$SD_Dep_Sham_Baseline*df.dep$SD_Dep_Sham_End) # Meta-analysis for depression changes rma_dep4 <- metacont(n.e=N_Active, mean.e = m.dep.ch.act ,sd.e = sd.dep.ch.act , n.c = N_Sham, mean.c = m.dep.ch.shm, sd.c = sd.dep.ch.shm, studlab = StudyID, byvar = Target, sm = "SMD", method.smd= "Hedges", data=df.dep, comb.fixed = FALSE) rma_dep4 meta::forest(rma_dep4, label.right = "Favors Active", label.left= "Favors Sham", digits = 3, digits.se=3, lab.e="Active TMS", lab.c="Sham TMS", hetstat= FALSE) df.dep$ES_Dep <- rma_dep4$TE df.dep$V_Dep <- (rma_dep4$seTE)^2 # Calculating anxiety mean changes c <- 0.5 # Keeping the correlation between the pre-post scores to be this df.anx$m.anx.ch.act <- df.anx$Mean_Anx_Active_Baseline - df.anx$Mean_Anx_Active_End df.anx$sd.anx.ch.act <- sqrt(df.anx$SD_Anx_Active_Baseline^2 + df.anx$SD_Anx_Active_End^2 + 2*c*df.anx$SD_Anx_Active_Baseline*df.anx$SD_Anx_Active_End) df.anx$m.anx.ch.shm <- df.anx$Mean_Anx_Sham_Baseline - df.anx$Mean_Anx_Sham_End df.anx$sd.anx.ch.shm <-sqrt(df.anx$SD_Anx_Sham_Baseline^2 + df.anx$SD_Anx_Sham_End^2 + 2*c*df.anx$SD_Anx_Sham_Baseline*df.anx$SD_Anx_Sham_End) # Meta-analysis for anxiety changes rma_anx4 <- metacont(n.e=N_Active, mean.e = m.anx.ch.act ,sd.e = sd.anx.ch.act , n.c = N_Sham, mean.c = m.anx.ch.shm, sd.c = sd.anx.ch.shm, studlab = StudyID, byvar = Target, sm = "SMD", method.smd= "Hedges", data=df.anx, comb.fixed = FALSE) rma_anx4 meta::forest(rma_anx4, label.right = "Favors Active", hetstat= FALSE, label.left= "Favors Sham", digits = 3, digits.se=3) df.anx$ES_Anx <- rma_anx4$TE df.anx$V_Anx <- (rma_anx4$seTE)^2 # Calculating YBOCS mean changes c <- 0.5 # Keeping the correlation between the pre-post scores to be this df.YBOCS$m.YBOCS.ch.act <- df.YBOCS$Mean_YBOCS_Active_Baseline - df.YBOCS$Mean_YBOCS_Active_End df.YBOCS$sd.YBOCS.ch.act <- sqrt(df.YBOCS$SD_YBOCS_Active_Baseline^2 + df.YBOCS$SD_YBOCS_Active_End^2 + 2*c*df.YBOCS$SD_YBOCS_Active_Baseline*df.YBOCS$SD_YBOCS_Active_End) df.YBOCS$m.YBOCS.ch.shm <- df.YBOCS$Mean_YBOCS_Sham_Baseline - df.YBOCS$Mean_YBOCS_Sham_End df.YBOCS$sd.YBOCS.ch.shm <-sqrt(df.YBOCS$SD_YBOCS_Sham_Baseline^2 + df.YBOCS$SD_YBOCS_Sham_End^2 + 2*c*df.YBOCS$SD_YBOCS_Sham_Baseline*df.YBOCS$SD_YBOCS_Sham_End) # Meta-analysis for YBOCS changes rma_YBOCS4 <- metacont(n.e=N_Active, mean.e = m.YBOCS.ch.act ,sd.e = sd.YBOCS.ch.act , n.c = N_Sham, mean.c = m.YBOCS.ch.shm, sd.c = sd.YBOCS.ch.shm, studlab = StudyID, byvar = Target, sm = "SMD", method.smd= "Hedges", data=df.YBOCS, comb.fixed = FALSE, tau.common = FALSE) rma_YBOCS4 meta::forest(rma_YBOCS4, label.right = "Favors Active", label.left= "Favors Sham", leftlabs = c("Study", "N", "Mean Change", "SD Change", "N", "Mean Change", "SD Change"), digits = 3, digits.se=3, hetstat= FALSE) df.YBOCS$ES_YBOCS <- rma_YBOCS4$TE df.YBOCS$V_YBOCS <- (rma_YBOCS4$seTE)^2 # Meta-Regression for anxiety predicting YBOCS change df.YBOCS.anx <- merge(df.YBOCS, select(df.anx, StudyID, contains("anx.ch"), contains("_Anx")), by="StudyID") df.YBOCS.anx$es.ch.anx <- df.YBOCS.anx$es.anx.ch.shm - df.YBOCS.anx$es.anx.ch.act rma_YBOCS4 <- metacont(n.e=N_Active, mean.e = m.YBOCS.ch.act ,sd.e = sd.YBOCS.ch.act , n.c = N_Sham, mean.c = m.YBOCS.ch.shm, sd.c = sd.YBOCS.ch.shm, studlab = StudyID, byvar = Target, sm = "SMD", method.smd= "Hedges", data=df.YBOCS.anx, comb.fixed = FALSE, tau.common = FALSE) rma_YBOCS4.anxreg <- metareg(rma_YBOCS4, ~ ES_Anx) rma_YBOCS4.anxreg bubble(rma_YBOCS4.anxreg, studlab = TRUE, cex.studlab = 0.7, pos.studlab = 1, xlab = "Effect Size of Anxiety Change", ylab="Effect Size of YBOCS Change") # Meta-Regression for depression predicting YBOCS change df.YBOCS.dep <- merge(df.YBOCS, select(df.dep, StudyID, contains("dep.ch"),contains("_Dep")), by="StudyID") df.YBOCS.dep$es.ch.dep <- df.YBOCS.dep$es.dep.ch.shm - df.YBOCS.dep$es.dep.ch.act rma_YBOCS5 <- metacont(n.e=N_Active, mean.e = m.YBOCS.ch.act ,sd.e = sd.YBOCS.ch.act , n.c = N_Sham, mean.c = m.YBOCS.ch.shm, sd.c = sd.YBOCS.ch.shm, studlab = StudyID, byvar = Target, sm = "SMD", method.smd= "Hedges", data=df.YBOCS.dep, comb.fixed = FALSE, tau.common = FALSE) rma_YBOCS5.depreg <- metareg(rma_YBOCS5, ~ ES_Dep) rma_YBOCS5.depreg bubble(rma_YBOCS5.depreg, studlab = TRUE, cex.studlab = 0.7, pos.studlab = 1, xlab = "Effect Size of Depression Change", ylab="Effect Size of YBOCS Change") #MetaSEM library(metaSEM) detach("package:meta", unload = TRUE) df.YBOCS.anx$cov.YBOCS.anx <- with(df.YBOCS.anx, sqrt(V_Anx) * sqrt(V_YBOCS) * 0.8) anxSEM <- meta(y = cbind(ES_Anx, ES_YBOCS), v = cbind(V_Anx, cov.YBOCS.anx, V_YBOCS), data = df.YBOCS.anx) rerun(anxSEM) summary(anxSEM) df.YBOCS.anx[,c("ES_Anx", "V_Anx", "ES_YBOCS", "V_YBOCS", "cov.YBOCS.anx")] with(df.YBOCS.anx, cor.test(ES_Anx, ES_YBOCS)) # Transforming the data into correlation matrices cormat <- function(x, y) {cor.anx.ybocs= y es.ybocs =x[,"ES_YBOCS"] es.anxiety = x[,"ES_Anx"] tmp <- matrix(c(1,cor.anx.ybocs,es.ybocs, cor.anx.ybocs,1,es.anxiety, es.ybocs,es.anxiety,1), nrow=3, byrow=T, dimnames = list(c("YBOCS", "Anxiety", "Treatment"), c("YBOCS", "Anxiety", "Treatment"))) return(tmp)} corlist <- list() for(i in 1:nrow(df.YBOCS.anx)) {corlist$data[[i]] <- cormat(df.YBOCS.anx[i,], 0.6)} corlist$N <- df.YBOCS.anx$N_Active + df.YBOCS.anx$N_Sham which(is.pd(corlist$data)==FALSE) corlist2 <- lapply(corlist, function(x) x[-c(4)]) obs.vars1 <- c("YBOCS", "Anxiety", "Treatment") ### Stage 1 analysis ## Random-effects model random1 <- tssem1(corlist2$data, corlist2$N, method="REM") rerun(random1) summary(random1) ## Average correlation matrix under a random-effects model averageR <- vec2symMat(coef(random1, select="fixed"), diag = FALSE) dimnames(averageR) <- list(obs.vars1, obs.vars1) averageR ## Heterogeneity variances of the random-effects coef(random1, select="random") ## Proposed model in lavaan syntax model1 <- "Anxiety ~ c*Treatment + b*YBOCS YBOCS ~ a*Treatment Treatment ~~ 1*Treatment" plot(model1) ## Convert the lavaan syntax to RAM specification used in metaSEM RAM1 <- lavaan2RAM(model1, obs.variables=obs.vars1) RAM1 tssem.fit <- tssem2(random1, RAM=RAM1, intervals.type = "LB", mx.algebras = list(ind=mxAlgebra(a*b, name="ind"), dir=mxAlgebra(c, name="dir"))) summary(tssem.fit) plot(tssem.fit, color="green") # Depression df.YBOCS.dep$cov.YBOCS.dep <- with(df.YBOCS.dep, sqrt(V_Dep) * sqrt(V_YBOCS) * 0.8) df.YBOCS.dep[,c("ES_Dep", "V_Dep", "ES_YBOCS", "V_YBOCS", "cov.YBOCS.dep")] with(df.YBOCS.dep, cor.test(ES_Dep, ES_YBOCS)) # Transforming the data into correlation matrices cormat <- function(x, y) {cor.dep.ybocs= y es.ybocs =x[,"ES_YBOCS"] es.dep = x[,"ES_Dep"] tmp <- matrix(c(1,cor.dep.ybocs,es.ybocs, cor.dep.ybocs,1,es.dep, es.ybocs,es.dep,1), nrow=3, byrow=T, dimnames = list(c("YBOCS", "Depression", "Treatment"), c("YBOCS", "Depression", "Treatment"))) return(tmp)} corlist <- list() for(i in 1:nrow(df.YBOCS.dep)) {corlist$data[[i]] <- cormat(df.YBOCS.dep[i,], 0.2)} corlist$N <- df.YBOCS.dep$N_Active + df.YBOCS.dep$N_Sham which(is.pd(corlist$data)==FALSE) corlist2 <- lapply(corlist, function(x) x[-c(4,6)]) obs.vars1 <- c("YBOCS", "Depression", "Treatment") ### Stage 1 analysis ## Random-effects model random1 <- tssem1(corlist2$data, corlist2$N, method="REM") rerun(random1) summary(random1) ## Average correlation matrix under a random-effects model averageR <- vec2symMat(coef(random1, select="fixed"), diag = FALSE) dimnames(averageR) <- list(obs.vars1, obs.vars1) averageR ## Heterogeneity variances of the random-effects coef(random1, select="random") ## Proposed model in lavaan syntax model1 <- "Depression ~ c*Treatment + b*YBOCS YBOCS ~ a*Treatment Treatment ~~ 1*Treatment" plot(model1) ## Convert the lavaan syntax to RAM specification used in metaSEM RAM1 <- lavaan2RAM(model1, obs.variables=obs.vars1) RAM1 tssem.fit <- tssem2(random1, RAM=RAM1, intervals.type = "LB", mx.algebras = list(ind=mxAlgebra(a*b, name="ind"), dir=mxAlgebra(c, name="dir"))) summary(tssem.fit) plot(tssem.fit, color="green")