### Remove any old objects from the R environment rm(list=ls()) # ------------------------------------------------------------- #Access libraries and accessing files # ------------------------------------------------------------- ### Access libraries library(dplyr) library(tidyverse) library(meta) library(metafor) library(car) library(rempsyc) library(ggplot2) library(ggsignif) options(stringsAsFactors = FALSE) ###read data files library(readxl) dat3 <- read_excel("~/Documents/Meta-Analysis/Data.xlsx") #exchange with your path # ------------------------------------------------------------- # Variable names and conversion from string to numeric # ------------------------------------------------------------- # "ID" (study name) # "EN" (N in intervention group) # "EMean" (mean in intervention group) # "ESD" (SD in intervention group) # "EN" (N in control group) # "CMean" (mean of control group outcome) # "CSD" (SD of control group outcome) ### Convert variables into numeric dat3$ESD <- as.numeric(dat3$ESD) dat3$CSD <- as.numeric(dat3$CSD) dat3$EMean <- as.numeric(dat3$EMean) dat3$CMean <- as.numeric(dat3$CMean) dat3$EN <- as.numeric(dat3$EN) dat3$CN <- as.numeric(dat3$CN) dat3$t <- as.numeric(dat3$t) dat3$tDCS.intensity <- as.numeric(dat3$tDCS.intensity) dat3$electrode.size <- as.numeric(dat3$electrode.size) dat3$Pubyear <- as.numeric(dat3$Pubyear) dat3$Duration <- as.numeric(dat3$Duration) dat3$Tarlocation <- factor(dat3$Tarlocation) dat3$task <- factor(dat3$task) dat3$tDCS.polarity <- factor(dat3$tDCS.polarity) dat3$return <- factor(dat3$return) dat3$timing <- factor(dat3$timing) dat3$ID <- factor(dat3$ID) # ------------------------------------------------------------- # Computing additional variables # ------------------------------------------------------------- #current density dat3 <- dat3 %>% mutate(current.density_calc = dat3$tDCS.intensity / dat3$electrode.size) dat3$current.density_calc <- as.numeric(dat3$current.density_calc) #####Descriptives####### Summarytable <-cbind(dat3$ID, dat3$EN, dat3$CN, dat3$control.condition, dat3$Tarlocation, dat3$return, dat3$current.density_calc, dat3$Duration, dat3$timing, dat3$task) #ADD OUTCOME AND FINDINGS # ------------------------------------------------------------- # Pre-processing: Imputation of correlation coefficient for crossover trials # ------------------------------------------------------------- # ------------------------------------------------------------- # Dealing with missing SD differences --> Calculating the standard deviation difference by using the t-statistic # ------------------------------------------------------------- # Function to calculate sd only if all required values are available dat3 <- dat3 %>% mutate (sd_diff = ifelse(!is.na(t), (dat3$EMean - dat3$CMean) / ((dat3$t) / sqrt(dat3$EN)), NA)) # # calculating the correlation coefficient for each observation in the filtered data frame # dat3_filtered <- dat3 %>% filter(!is.na(sd_diff) & EN > 0) if (nrow(dat3_filtered) > 0) {} dat3_filtered <- dat3_filtered %>% mutate (CORRCOE = (ESD^2 + CSD^2 - sd_diff^2) / (2 * ESD * CSD)) #Formula from Cochrane Mean_CORRCOE <- mean(dat3_filtered$CORRCOE) ### calculating imputed CORRCOE from Dormal study, optional #imp <- dat3 %>% filter(!is.na(DiffESD)) #imp <- imp %>% mutate (CORRCOE1 = (ESD^2 + CSD^2 - DiffESD^2) / (2 * ESD * CSD)) #The difference SD is based on the summary data in the article. Calculating the correlation coefficient based on the value can be helpful to validate the method using the t-formula to impute missing difference SD # ------------------------------------------------------------- # Descriptives: # ------------------------------------------------------------- Total_sample_size <- sum(dat3$EN) + sum(dat3$CN) # ------------------------------------------------------------- # Effect size calculation: # ------------------------------------------------------------- ###Dividing sample size of the sham by number of trial arms to account for multiple comparison dat3$CN <- (dat3$CN/dat3$comp) # splitting data file in parallel and crossover studies for calculating the effect sizes independently with different formulas datP <- dat3 %>% filter(control.condition == "between-subjects") datC <- dat3 %>% filter(control.condition == "within") # ------------------------------------------------------------- # Effect size calculation: for parallel trials: # ------------------------------------------------------------- ###lnCVR calculation for parallel trials including correction for different sample sizes from Senior, A. M., Viechtbauer, W., & Nakagawa, S. (2020). Revisiting and expanding the meta‐analysis of variation: The log coefficient of variation ratio. Research Synthesis Methods, 11(4), 553-567. ####Point estimator datP$CVRyi <- log((datP$ESD / datP$EMean) / (datP$CSD / datP$CMean)) datP$CVRyi <- datP$CVRyi + (1 / (2 * (datP$EN - 1))) - (1 / (2 * (datP$CN - 1))) datP$CVRyi <- datP$CVRyi + 0.5 * (((datP$CSD^2) / ((datP$CN * datP$CMean^2))) - ((datP$ESD^2) / ((datP$EN * datP$EMean^2)))) ####Sampling distribution datP$CVRvi <- datP$CSD^2 / (datP$CN * datP$CMean^2) + datP$CSD^4 / (2 * datP$CN^2 * datP$CMean^4) + 1 / (2 * (datP$CN - 1)) + 1 / (2 * (datP$CN - 1)^2) + datP$ESD^2 / (datP$EN * datP$EMean^2) + datP$ESD^4 / (2 * datP$EN^2 * datP$EMean^4) + 1 / (2 * (datP$EN - 1)) + 1 / (2 * (datP$EN - 1)^2) # ------------------------------------------------------------- # Effect size calculation: Crossover studies: with CORRCOE # ------------------------------------------------------------- # CORRCOE calculation Crossover trials adapted from Senior, A. M., Viechtbauer, W., & Nakagawa, S. (2020). Revisiting and expanding the meta‐analysis of variation: The log coefficient of variation ratio. Research Synthesis Methods, 11(4), 553-567. ### creating N variable for formula below (N = EN = CN) in cross over trails and determining correlation between active and sham ###### datC <- datC %>% mutate(N = (datC$EN + datC$CN)/2) ##Manually calculating corrected lnCVR (yic) and corrected sampling distribution (vic) datC$CVRyi <- log((datC$ESD / datC$EMean) / (datC$CSD / datC$CMean)) + 0.5 * (datC$CSD^2 / (datC$CN * datC$CMean^2) - datC$ESD^2 / (datC$EN * datC$EMean^2)) datC$CVRvi <- datC$CSD^2 / (datC$N * datC$CMean^2) + datC$ESD^2 / (datC$N * datC$EMean^2) - Mean_CORRCOE * ((2 * datC$CSD * datC$ESD) / (datC$N * datC$CMean * datC$EMean)) + 1 / (datC$N - 1) - Mean_CORRCOE^2 * (1 / (datC$N - 1)) # ------------------------------------------------------------- # For sensitivity analysis # ------------------------------------------------------------- # ------------------------------------------------------------- # Effect size calculation: Crossover studies: with 0.3 # ------------------------------------------------------------- #crossover datC$CVRyi0.3 <- log((datC$ESD / datC$EMean) / (datC$CSD / datC$CMean)) + 0.5 * (datC$CSD^2 / (datC$CN * datC$CMean^2) - datC$ESD^2 / (datC$EN * datC$EMean^2)) datC$CVRvi0.3 <- datC$CSD^2 / (datC$N * datC$CMean^2) + datC$ESD^2 / (datC$N * datC$EMean^2) - 0.3 * ((2 * datC$CSD * datC$ESD) / (datC$N * datC$CMean * datC$EMean)) + 1 / (datC$N - 1) - 0.3^2 * (1 / (datC$N - 1)) #Parallel datP$CVRyi0.3 <- log((datP$ESD / datP$EMean) / (datP$CSD / datP$CMean)) datP$CVRyi0.3 <- datP$CVRyi + (1 / (2 * (datP$EN - 1))) - (1 / (2 * (datP$CN - 1))) datP$CVRyi0.3 <- datP$CVRyi + 0.5 * (((datP$CSD^2) / ((datP$CN * datP$CMean^2))) - ((datP$ESD^2) / ((datP$EN * datP$EMean^2)))) datP$CVRvi0.3 <- datP$CSD^2 / (datP$CN * datP$CMean^2) + datP$CSD^4 / (2 * datP$CN^2 * datP$CMean^4) + 1 / (2 * (datP$CN - 1)) + 1 / (2 * (datP$CN - 1)^2) + datP$ESD^2 / (datP$EN * datP$EMean^2) + datP$ESD^4 / (2 * datP$EN^2 * datP$EMean^4) + 1 / (2 * (datP$EN - 1)) + 1 / (2 * (datP$EN - 1)^2) # ------------------------------------------------------------- # Effect size calculation: Crossover studies: with 0.8 # ------------------------------------------------------------- #crossover datC$CVRyi0.8 <- log((datC$ESD / datC$EMean) / (datC$CSD / datC$CMean)) + 0.5 * (datC$CSD^2 / (datC$CN * datC$CMean^2) - datC$ESD^2 / (datC$EN * datC$EMean^2)) datC$CVRvi0.8 <- datC$CSD^2 / (datC$N * datC$CMean^2) + datC$ESD^2 / (datC$N * datC$EMean^2) - 0.8 * ((2 * datC$CSD * datC$ESD) / (datC$N * datC$CMean * datC$EMean)) + 1 / (datC$N - 1) - 0.8^2 * (1 / (datC$N - 1)) #Parallel datP$CVRyi0.8 <- log((datP$ESD / datP$EMean) / (datP$CSD / datP$CMean)) datP$CVRyi0.8 <- datP$CVRyi + (1 / (2 * (datP$EN - 1))) - (1 / (2 * (datP$CN - 1))) datP$CVRyi0.8 <- datP$CVRyi + 0.5 * (((datP$CSD^2) / ((datP$CN * datP$CMean^2))) - ((datP$ESD^2) / ((datP$EN * datP$EMean^2)))) datP$CVRvi0.8 <- datP$CSD^2 / (datP$CN * datP$CMean^2) + datP$CSD^4 / (2 * datP$CN^2 * datP$CMean^4) + 1 / (2 * (datP$CN - 1)) + 1 / (2 * (datP$CN - 1)^2) + datP$ESD^2 / (datP$EN * datP$EMean^2) + datP$ESD^4 / (2 * datP$EN^2 * datP$EMean^4) + 1 / (2 * (datP$EN - 1)) + 1 / (2 * (datP$EN - 1)^2) # ------------------------------------------------------------- # ------------------------------------------------------------- # ------------------------------------------------------------- # ------------------------------------------------------------- # MAIN ANALYSIS: # ------------------------------------------------------------- dat <- bind_rows(datP, datC) dat <- escalc(measure = "SMD", m1i=EMean, sd1i=ESD, n1i=EN, m2i=CMean, sd2i=CSD, n2i=CN, data=dat, var.names=c("SMD_yi","SMD_vi"), append = TRUE) ############# Caluculating the number of studies that had a lnCVR < 1########### # Overlay Gaussian curve mu <- mean(dat$CVRyi) sigma <- sd(dat$CVRyi) curve(dnorm(x, mean = mu, sd = sigma), col = "red", lwd = 2, add = TRUE) # Plot histogram of CVRyi hist(dat$CVRyi, freq = FALSE, main = "Histogram of CVRyi with Gaussian Curve", xlab = "CVRyi", ylab = "Density", col = "lightblue") # Overlay Gaussian curve mu <- mean(dat$CVRyi) sigma <- sd(dat$CVRyi) curve(dnorm(x, mean = mu, sd = sigma), col = "red", lwd = 2, add = TRUE) plot(y = dat$CVRvi, x = dat$CVRyi) qqnorm(dat$CVRyi, main = "QQ Plot of CVRyi against Normal Distribution") qqline(dat$CVRyi, col = "red") ### Main families 1 and 2 MLMA.1 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1 | control.condition, ~ 1 | Pubyear), data=dat, method="REML") summary(MLMA.1) funnel(MLMA.1, xlab = "lnCVR") results.Main <- coef(summary(MLMA.1)) exp(results.Main$estimate) exp(results.Main$ci.lb) exp(results.Main$ci.ub) Main_forest <- forest(MLMA.1, atransf = exp, at = log(c(0.1, 1, 5, 10)), addpred = TRUE, header = TRUE, slab = ID, xlab = "CVR", xlim = c(-5, 5), alim = c(-7, 4), ylim = c(91.2785, 0.1), cex.axis = 1, mlab = "", #ilab = cbind(active_N, sham_N), ilab.xpos = c(-4.5, -3.5)) op <- par(cex = 0.75, font = 2) text(c(-4.5, -3.5), Main_forest$k + 2, c("Active", "Sham"), cex = 1, font = 1) par(op) datexp <- dat %>% mutate(CVRyi = exp(CVRyi)) datPerc <- datexp %>% mutate(datcutoff = ifelse(CVRyi < 1, "Below 1", "Above 1")) Below_1 <- sum(datPerc$datcutoff == "Below 1") Above_1 <- sum(datPerc$datcutoff == "Above 1") regtest(SMD_yi, SMD_vi, data=dat) ranktest(SMD_yi, SMD_vi, data=dat) funnel(MLMA.1, xlab = "lnCVR") # Assuming 'rma.mv_object' is your multivariate random-effects meta-analysis object effects <- MLMA.1$b se <- sqrt(diag(MLMA.1$V)) corr <- MLMA.1$rho prec <- sum(1/se) pet_peese_test <- lm(effects ~ prec + sqrt(1 - corr^2)/prec) summary(pet_peese_test) # Perform PET-PEESE pet_peese_test <- lm(effects ~ 1/prec) # Print the results summary(pet_peese_test) regtest(SMD_yi, SMD_vi, data=dat) ranktest(SMD_yi, SMD_vi, data=dat) funnel(taf, legend=TRUE) taf regtest(main) ################################################################################ # stepwise procedure for CVR to include seperate moderators, based on Schroeder et al. 2020 MLMA.SW1 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(task), data=dat, method="REML") summary(MLMA.SW1) MLMA.SW2 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(tDCS.polarity), data=dat, method="REML") summary(MLMA.SW2) MLMA.SW3 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tarlocation), data=dat, method="REML") summary(MLMA.SW3) MLMA.SW3.2 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ factor(Tlocation), data=dat, method="REML") summary(MLMA.SW3.2) MLMA.SW4 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ factor(return), data=dat, method="REML") summary(MLMA.SW4) MLMA.SW5 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ factor(timing), data=dat, method="REML") summary(MLMA.SW5) MLMA.SW6 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ current.density_calc, data=dat, method="REML") summary(MLMA.SW6) MLMA.SW7 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ Duration, data=dat, method="REML") summary(MLMA.SW7) MLMA.SW8 <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear) , mods = ~ factor(blinding), data=dat, method="REML") summary(MLMA.SW8) MLMA.SWX <- rma.mv(dat$CVRyi, dat$CVRvi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tlocation) * Duration, data=dat, method="REML") summary(MLMA.SWX) ### --> does not work # stepwise procedure for SMD to include seperate moderators, based on Schroeder et al. 2020 MLMA.SW.SMD <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), data=dat, method="REML") summary(MLMA.SW.SMD) funnel(MLMA.SW.SMD, xlab = "Hedges’ g", xlim = c(-3.25, 2)) ranktest(MLMA.SW.SMD) regtest(dat$SMD_yi, dat$SMD_vi, data = dat) Gregg <- rma(dat$SMD_yi, dat$SMD_vi, random = ~ 1 | ID / control.condition / Pubyear, data=dat, method="REML") summary(Gregg) trimfill(Gregg) SMD_forest <- forest(MLMA.SW.SMD, at = log(c(0.1, 1, 5, 10)), addpred = TRUE, header = TRUE, slab = ID, xlab = "hedges g", xlim = c(-6, 5), alim = c(-8, 5), ylim = c(91.2785, 0.1), cex.axis = 1, mlab = "Overall hedges g", ilab.xpos = c(-4.5, -3.5)) MLMA.SW.SMD.1 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(task) -1, data=dat, method="REML") summary(MLMA.SW.SMD.1) MLMA.SW.SMD.2 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(tDCS.polarity) -1, data=dat, method="REML") summary(MLMA.SW.SMD.2) MLMA.SW.SMD.3 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tarlocation) -1 , data=dat, method="REML") summary(MLMA.SW.SMD.3) MLMA.SW.SMD.3.1 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tlocation) -1, data=dat, method="REML") summary(MLMA.SW.SMD.3.1) MLMA.SW.SMD.4 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(return) -1 , data=dat, method="REML") summary(MLMA.SW.SMD.4) MLMA.SW.SMD.4.2 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(returnSchroeder) -1 , data=dat, method="REML") summary(MLMA.SW.SMD.4.2) MLMA.SW.SMD.5 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(timing) -1, data=dat, method="REML") summary(MLMA.SW.SMD.5) MLMA.SW.SMD.6 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ current.density_calc, data=dat, method="REML") summary(MLMA.SW.SMD.6) MLMA.SW.SMD.7 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ Duration, data=dat, method="REML") summary(MLMA.SW.SMD.7) MLMA.SW.SMD.8 <- rma.mv(dat$SMD_yi, dat$SMD_vi, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(blinding), data=dat, method="REML") summary(MLMA.SW.SMD.8) ###table library(apaTables) results.SMD1 <- coef(summary(MLMA.SW.SMD.1)) results.SMD2 <- coef(summary(MLMA.SW.SMD.2)) results.SMD3 <- coef(summary(MLMA.SW.SMD.3)) results.SMD4 <- coef(summary(MLMA.SW.SMD.4)) results.SMD5 <- coef(summary(MLMA.SW.SMD.5[, "estimate"])) results.SMD6 <- coef(summary(MLMA.SW.SMD.6[, "estimate"])) results.SMD7 <- coef(summary(MLMA.SW.SMD.7[, "estimate"])) results.SMD8 <- coef(summary(MLMA.SW.SMD.8[, "estimate"])) results.SMD9 <- coef(summary(MLMA.SW.SMD.9[, "estimate"])) Level <- nrowbind("Contralateral", "HD", "Ipsilateral", "supraorbital", "extracranial") k <- nrow(results.SMD4) # Number of coefficients p_value <- results.SMD4[, "estimate"] # p-values z_value <- results.SMD4[, "zval"] # z-values estimate <- results.SMD4[, "pval"] reference.table <- cbind(Level, k, estimate, z_value, p_value) reference.table <- as.data.frame(reference.table) nice_table(reference.table) healthy <- rma.mv(CVRyi, CVRvi, subset = (population == "healthy"), data = dat ) healthy_forest <- forest(healthy, atransf = exp, at = log(c(0.1, 1, 5, 10)), addpred = TRUE, header = TRUE, slab = ID, xlab = "CVR", xlim = c(-7, 5), alim = c(-8, 5), ylim = c(91.2785, 0.1), cex.axis = 1, mlab = "", ilab = cbind(active_N, sham_N), ilab.xpos = c(-4.5, -3.5)) ##ggplot for SMD task results <- coef(summary(MLMA.SW.SMD.1)) # Create a data frame for plotting plot_task <- data.frame( Group = c("GNGT", "SST"), Effect_Size = results[, "estimate"], P = results[, "pval"], CI_Lower = results[, "ci.lb"], CI_Upper = results[, "ci.ub"] ) # Create the ggplot ggplot(plot_task, aes(x = Group, y = Effect_Size)) + geom_bar(stat = "identity", aes(fill = Group)) + geom_errorbar(aes(ymax = CI_Upper, ymin = CI_Lower), width = 0.25, position = position_dodge(width = 0.1), size = 0.25) + labs( x = "Task", y = "Hedges’ g") + theme_bw() + theme(legend.position = "none") + theme( text = element_text(family = "Times New Roman"), plot.title = element_text(size = 12, hjust = 0.5), axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), legend.title = element_text(size = 10), legend.text = element_text(size = 10), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), ) + coord_cartesian(ylim = c(-0.5, 0.5)) results123 <- coef(summary(SSTMV)) # Create a data frame for plotting plot_task2 <- data.frame( Group = c("SST"), Effect_Size = results123[, "estimate"], P = results123[, "pval"], CI_Lower = results123[, "ci.lb"], CI_Upper = results123[, "ci.ub"] ) # Create the ggplot ggplot(plot_task2, aes(x = Group, y = Effect_Size)) + geom_bar(stat = "identity", aes(fill = Group)) + geom_errorbar(aes(ymax = CI_Upper, ymin = CI_Lower), width = 0.25, position = position_dodge(width = 0.1), size = 0.25) + labs(title = "Mean difference between SST and GNGT", x = "SST", y = "headges g") + theme_bw() + theme(legend.position = "none") + theme( text = element_text(family = "Times New Roman"), plot.title = element_text(size = 12, hjust = 0.5), axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), legend.title = element_text(size = 10), legend.text = element_text(size = 10), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), ) + coord_cartesian(ylim = c(-0.5, 0.5)) + # Set y-axis limits geom_signif(comparisons = list(c("SST")), map_signif_level = TRUE) # ------------------------------------------------------------- # SUB ANALYSIS Population : # ------------------------------------------------------------- dat$population <- fct_collapse(dat$population, "healthy" = c("healthy"), "patient" = c("MDD", "veterans", "alcohol", "tabacco", "healthy and PWS", "other patients", "Parkinson")) # ------------------------------------------------------------- # SENSITIVITY ANALYSIS: # ------------------------------------------------------------- # ------------------------------------------------------------- # Effect size calculation: Crossover studies: corr = 0.3 # ------------------------------------------------------------- # Crossover trials with 0.3 MLMA.0.3 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), data=dat, method="REML",) summary(MLMA.0.3) MLMA.SW1.03 <- (rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(task), data=dat, method="REML")) summary(MLMA.SW1.03) MLMA.SW2.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(tDCS.polarity), data=dat, method="REML") summary(MLMA.SW2.03) MLMA.SW3.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tlocation), data=dat, method="REML") summary(MLMA.SW3.03) MLMA.SW4.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(return), data=dat, method="REML") summary(MLMA.SW4.03) MLMA.SW5.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(timing), data=dat, method="REML") summary(MLMA.SW5.03) MLMA.SW6.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ current.density_calc, data=dat, method="REML") summary(MLMA.SW6.03) MLMA.SW7.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ Duration, data=dat, method="REML") summary(MLMA.SW7.03) MLMA.SW8.03 <- rma.mv(dat$CVRyi0.3, dat$CVRvi0.3, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(blinding), data=dat, method="REML") summary(MLMA.SW8.03) # Crossover trials with 0.8 MLMA.0.8 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = ~ 1 | ID, data=dat, method="REML",) summary(MLMA.0.8) MLMA.SW1.08 <- (rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(task), data=dat, method="REML")) summary(MLMA.SW1.08) MLMA.SW2.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(tDCS.polarity), data=dat, method="REML") summary(MLMA.SW2.08) MLMA.SW3.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(Tlocation), data=dat, method="REML") summary(MLMA.SW3.08) MLMA.SW4.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(return), data=dat, method="REML") summary(MLMA.SW4.08) MLMA.SW5.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(timing), data=dat, method="REML") summary(MLMA.SW5.08) MLMA.SW6.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ current.density_calc, data=dat, method="REML") summary(MLMA.SW6.08) MLMA.SW7.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ Duration, data=dat, method="REML") summary(MLMA.SW7.08) MLMA.SW8.08 <- rma.mv(dat$CVRyi0.8, dat$CVRvi0.8, random = list (~ 1 | ID, ~ 1| control.condition, ~ 1 | Pubyear), mods = ~ factor(blinding), data=dat, method="REML") summary(MLMA.SW8.08) datBias <- dat %>% filter(Bias == "Low" | Bias == "Moderate") MLMA.BiasCVR <- rma.mv(datBias$CVRyi, datBias$CVRvi, random = ~ 1 | ID, data=datBias, method="REML",) summary(MLMA.BiasCVR) MLMA.BiasSMD <- rma.mv(datBias$SMD_yi, datBias$SMD_vi, random = ~ 1 | ID, data=datBias, method="REML",) summary(MLMA.BiasSMD) new.df <- dat[dat$ID != "Chen (2021)", ] ochen <- rma.mv(new.df$SMD_yi, new.df$SMD_vi, random = list (~ 1 | ID, ~ 1 | control.condition, ~ 1 | Pubyear), mods = ~ factor(task) -1 , data=new.df, method="REML") summary(ochen) ############# Caluculating the number of studies that had a lnCVR < 1########### N_VR_below1 <- vrdat_trans %>% group_by(id) %>% mutate(av_VR = mean(zi)) %>% mutate(below_over_1 = ifelse(av_VR < 1, "below", "equal to 1 or over")) %>% group_by(below_over_1) %>% count(n_distinct(id)) #Transforming the values back from the log scale. From Homan et al. 2021 vrdat_trans <- summary(vrdat, trans = exp, digits = 2) #accounting for effect from the same study (size dependencies) escalc(measure=CVR, m1i = Final$Int_Mean, n1i = Final$Int_N, sd1i = Final$Int_SD, m2i = Finalt$Con_Mean, n2i = Final$Con_N, sd2i = Final$Con_SD ) mes_ss <- robu(formula =yi, data = data_mes, studynum =id, var.eff.size =vi, modelweights = "HIER", small =TRUE) #filter for size of CVR datAbove1 <- datPerc %>% filter(datcutoff == "Above 1") datBelow1 <- datPerc %>% filter(datcutoff == "Below 1") datAbove1 <- datAbove1 %>% mutate(standard_error_active_above1 =(datAbove1$ESD)/(sqrt(EN))) datBelow1 <- datBelow1 %>% mutate(standard_error_active_below1 =(datBelow1$ESD)/(sqrt(EN))) mean(datAbove1$standard_error_active_above1) mean(datBelow1$standard_error_active_below1) t.test(dat3$standard_error_sham, datAbove1$standard_error_active_above1) t.test(dat3$standard_error_sham, datBelow1$standard_error_active_below1) dat$EMean <- as.numeric(dat$EMean) dat$CMean <- as.numeric(dat$CMean) dat$ESD <- as.numeric(dat$ESD) dat$CSD <- as.numeric(dat$CSD) cor(dat$EMean, dat$ESD) plot(dat$EMean, dat$ESD) cor(dat$CMean, dat$CSD) plot(dat$CMean, dat$CSD)