Neuromodulation Meta-Analysis in Schizophrenia: Efficacy Analysis

1. Loading Required Packages

# Install packages if not already installed
packages <- c("metafor", "dplyr", "tidyr", "ggplot2", "readxl", "knitr", "kableExtra", "RColorBrewer")
for(pkg in packages) {
  if(!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
    library(pkg, character.only = TRUE)
  }
}

2. Importing Data

# Import the three sheets
baseline <- read_excel("Downloads/final_meta_analysis_data.xlsx", 
    sheet = "Baseline")
outcomes <- read_excel("Downloads/final_meta_analysis_data.xlsx", 
    sheet = "Outcomes")
#codes <- read_excel("Downloads/final_meta_analysis_data.xlsx", 
#    sheet = "Codes")

# Look at the structure of the imported data
str(baseline)
## tibble [256 × 16] (S3: tbl_df/tbl/data.frame)
##  $ Study            : chr [1:256] "Bais, 2014 - L" "Bais, 2014 - L" "Bais, 2014 - Bi" "Bais, 2014 - Bi" ...
##  $ Group            : chr [1:256] "Active" "Sham" "Active" "Sham" ...
##  $ N                : num [1:256] 16 16 15 16 12 12 13 12 13 12 ...
##  $ Age_Mean         : chr [1:256] "37.200000000000003" "37.299999999999997" "33.9" "37.299999999999997" ...
##  $ Age_SD           : chr [1:256] "14.9" "11.6" "9.1999999999999993" "11.6" ...
##  $ percentage_Male  : chr [1:256] "56.25" "62.5" "53.33" "62.5" ...
##  $ Target           : chr [1:256] "L-TPJ" "L-TPJ" "Bi-TPJ" "L-TPJ" ...
##  $ Lateralisation   : chr [1:256] "Uni-L" "Uni-L" "Bi-Seq" "Uni-L" ...
##  $ Targeting_Method : chr [1:256] "10-20EEG" "10-20EEG" "10-20EEG" "10-20EEG" ...
##  $ Technique        : chr [1:256] "LF-RTMS" "LF-RTMS" "LF-RTMS" "LF-RTMS" ...
##  $ Total_Time       : num [1:256] 240 240 240 240 25 25 500 500 500 500 ...
##  $ Total_Sessions   : num [1:256] 12 12 12 12 1 1 20 20 20 20 ...
##  $ Session_Frequency: chr [1:256] "high" "high" "high" "high" ...
##  $ Sham_Condition   : chr [1:256] "coil" "coil" "coil" "coil" ...
##  $ T0               : num [1:256] 0 0 0 0 0 0 0 0 0 0 ...
##  $ T1               : num [1:256] 8 8 8 8 0.0174 ...
str(outcomes)
## tibble [2,915 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Study          : chr [1:2915] "Bais, 2014 - L" "Bais, 2014 - L" "Bais, 2014 - L" "Bais, 2014 - L" ...
##  $ Group          : chr [1:2915] "Active" "Active" "Active" "Active" ...
##  $ Outcome_Measure: chr [1:2915] "P_PANSS" "P_PANSS" "N_PANSS" "N_PANSS" ...
##  $ Time_Point     : chr [1:2915] "T0" "T1" "T0" "T1" ...
##  $ Days           : num [1:2915] 0 8 0 8 0 8 0 8 0 8 ...
##  $ N              : num [1:2915] 16 16 16 16 16 16 16 16 16 14 ...
##  $ Mean           : num [1:2915] 16.3 15.1 15.1 14.5 30.1 ...
##  $ SD             : num [1:2915] 4.76 5.64 4.7 4.4 8.85 9.04 5.67 5.55 10.6 9.67 ...

3. Data Preprocessing

# Convert to lowercase for consistency
baseline$Study <- tolower(baseline$Study)
baseline$Group <- tolower(baseline$Group)
outcomes$Study <- tolower(outcomes$Study)
outcomes$Group <- tolower(outcomes$Group)

# Create unique identifiers for each study-group combination
baseline$study_group_id <- paste(baseline$Study, baseline$Group, sep = "_")
outcomes$study_group_id <- paste(outcomes$Study, outcomes$Group, sep = "_")

# Categorize outcome measures
outcomes <- outcomes %>%
  mutate(domain = case_when(
    grepl("^P_", Outcome_Measure) ~ "positive_symptoms",
    grepl("^N_", Outcome_Measure) ~ "negative_symptoms",
    grepl("^T_", Outcome_Measure) ~ "total_psychopathology",
    grepl("^G_", Outcome_Measure) ~ "general_psychopathology",
    grepl("^PS_", Outcome_Measure) ~ "processing_speed",
    grepl("^AV_", Outcome_Measure) ~ "attention_vigilance",
    grepl("^WM_", Outcome_Measure) ~ "working_memory",
    grepl("^VerL_", Outcome_Measure) ~ "verbal_learning",
    grepl("^VisL_", Outcome_Measure) ~ "visual_learning",
    grepl("^RPS_", Outcome_Measure) ~ "reasoning_problem_solving",
    grepl("^SC_", Outcome_Measure) ~ "social_cognition",
    grepl("^GC_", Outcome_Measure) ~ "global_cognition",
    TRUE ~ "other"
  ))

# Check the distribution of outcome domains
table(outcomes$domain)
## 
##       attention_vigilance   general_psychopathology          global_cognition 
##                       106                       216                       108 
##         negative_symptoms                     other         positive_symptoms 
##                       446                        91                       488 
##          processing_speed reasoning_problem_solving          social_cognition 
##                       188                       294                        78 
##     total_psychopathology           verbal_learning           visual_learning 
##                       304                       160                       142 
##            working_memory 
##                       294

4. Calculation of Effect Sizes

# Create a dataframe to store effect sizes
effect_sizes <- data.frame()

# Get unique study-group combinations
study_groups <- unique(outcomes$study_group_id)

# Identify active and sham groups
for(study in unique(outcomes$Study)) {
  # Get all groups for this study
  groups <- unique(outcomes$Group[outcomes$Study == study])
  
  # Identify which is active and which is sham
  active_groups <- groups[grepl("active", groups, ignore.case = TRUE)]
  sham_groups <- groups[grepl("sham", groups, ignore.case = TRUE)]
  
  for(active in active_groups) {
    for(sham in sham_groups) {
      # For each outcome measure
      for(outcome in unique(outcomes$Outcome_Measure)) {
        # Check if both groups have this outcome
        active_data <- outcomes %>% 
          filter(Study == study, Group == active, Outcome_Measure == outcome)
        
        sham_data <- outcomes %>% 
          filter(Study == study, Group == sham, Outcome_Measure == outcome)
        
        if(nrow(active_data) >= 2 && nrow(sham_data) >= 2) {
          # Get baseline and endpoint data
          active_t0 <- active_data %>% filter(Time_Point == "T0")
          active_t1 <- active_data %>% filter(Time_Point == "T1" | grepl("T1_Change", Time_Point))
          sham_t0 <- sham_data %>% filter(Time_Point == "T0")
          sham_t1 <- sham_data %>% filter(Time_Point == "T1" | grepl("T1_Change", Time_Point))
          
          # Only proceed if we have complete data
          if(nrow(active_t0) == 1 && nrow(active_t1) == 1 && 
             nrow(sham_t0) == 1 && nrow(sham_t1) == 1) {
            
            # Check if this is a change score
            is_change_score <- grepl("T1_Change", active_t1$Time_Point)
            
            # Calculate effect sizes based on whether it's a change score
            if(is_change_score) {
              # For change scores
              n1 <- as.numeric(active_t1$N)
              n2 <- as.numeric(sham_t1$N)
              m1 <- as.numeric(active_t1$Mean)
              m2 <- as.numeric(sham_t1$Mean)
              sd1 <- as.numeric(active_t1$SD)
              sd2 <- as.numeric(sham_t1$SD)
              
              # Skip if any data is missing
              if(any(is.na(c(n1, n2, m1, m2, sd1, sd2)))) next
              
              # Calculate Hedges' g
              es <- escalc(measure = "SMD", 
                         m1i = m1, m2i = m2, 
                         sd1i = sd1, sd2i = sd2, 
                         n1i = n1, n2i = n2)
            } else {
              # For pre-post data
              n1 <- as.numeric(active_t1$N)
              n2 <- as.numeric(sham_t1$N)
              m1_pre <- as.numeric(active_t0$Mean)
              m1_post <- as.numeric(active_t1$Mean)
              m2_pre <- as.numeric(sham_t0$Mean)
              m2_post <- as.numeric(sham_t1$Mean)
              sd1_pre <- as.numeric(active_t0$SD)
              sd1_post <- as.numeric(active_t1$SD)
              sd2_pre <- as.numeric(sham_t0$SD)
              sd2_post <- as.numeric(sham_t1$SD)
              
              # Skip if any data is missing
              if(any(is.na(c(n1, n2, m1_pre, m1_post, m2_pre, m2_post, 
                            sd1_pre, sd1_post, sd2_pre, sd2_post)))) next
              
              # Pre-post correlation (estimate 0.7 if not known)
              cor_pre_post <- 0.7
              
              # Calculate change scores and their SDs
              m1 <- m1_post - m1_pre
              m2 <- m2_post - m2_pre
              
              # Calculate SD of change scores
              sd1 <- sqrt(sd1_pre^2 + sd1_post^2 - 2*cor_pre_post*sd1_pre*sd1_post)
              sd2 <- sqrt(sd2_pre^2 + sd2_post^2 - 2*cor_pre_post*sd2_pre*sd2_post)
              
              # Calculate Hedges' g
              es <- escalc(measure = "SMD", 
                         m1i = m1, m2i = m2, 
                         sd1i = sd1, sd2i = sd2, 
                         n1i = n1, n2i = n2)
            }
            
            # Determine if lower scores are better (for direction adjustment)
            lower_is_better <- grepl("PANSS|SANS|SAPS|AHRS|BPRS|PSYRATS", outcome)
            
            # Adjust direction so positive effect size always means improvement
            if(lower_is_better) {
              es$yi <- -es$yi
            }
            
            # Get the domain of this outcome
            domain <- active_data$domain[1]
            
            # Get technique, target, and other metadata from baseline
            technique <- baseline$Technique[baseline$study_group_id == paste(study, active, sep="_")][1]
            target <- baseline$Target[baseline$study_group_id == paste(study, active, sep="_")][1]
            lateralisation <- baseline$Lateralisation[baseline$study_group_id == paste(study, active, sep="_")][1]
            targeting_method <- baseline$Targeting_Method[baseline$study_group_id == paste(study, active, sep="_")][1]
            total_time <- baseline$Total_Time[baseline$study_group_id == paste(study, active, sep="_")][1]
            total_sessions <- baseline$Total_Sessions[baseline$study_group_id == paste(study, active, sep="_")][1]
            session_frequency <- baseline$Session_Frequency[baseline$study_group_id == paste(study, active, sep="_")][1]
            
            # Add to effect sizes dataframe
            effect_sizes <- rbind(effect_sizes, data.frame(
              study = study,
              outcome = outcome,
              domain = domain,
              technique = technique,
              target = target,
              lateralisation = lateralisation,
              targeting_method = targeting_method,
              total_time = total_time,
              total_sessions = total_sessions,
              session_frequency = session_frequency,
              yi = es$yi,  # Effect size
              vi = es$vi,  # Variance
              n1 = n1,      # Sample size active
              n2 = n2,      # Sample size sham
              lower_is_better = lower_is_better,
              stringsAsFactors = FALSE
            ))
          }
        }
      }
    }
  }
}

# Check the resulting effect sizes
head(effect_sizes)
##            study outcome                  domain technique target
## 1 bais, 2014 - l P_PANSS       positive_symptoms   LF-RTMS  L-TPJ
## 2 bais, 2014 - l N_PANSS       negative_symptoms   LF-RTMS  L-TPJ
## 3 bais, 2014 - l G_PANSS general_psychopathology   LF-RTMS  L-TPJ
## 4 bais, 2014 - l  P_AHRS       positive_symptoms   LF-RTMS  L-TPJ
## 5 bais, 2014 - l P_PANAS       positive_symptoms   LF-RTMS  L-TPJ
## 6 bais, 2014 - l N_PANAS       negative_symptoms   LF-RTMS  L-TPJ
##   lateralisation targeting_method total_time total_sessions session_frequency
## 1          Uni-L         10-20EEG        240             12              high
## 2          Uni-L         10-20EEG        240             12              high
## 3          Uni-L         10-20EEG        240             12              high
## 4          Uni-L         10-20EEG        240             12              high
## 5          Uni-L         10-20EEG        240             12              high
## 6          Uni-L         10-20EEG        240             12              high
##           yi        vi n1 n2 lower_is_better
## 1  0.2912072 0.1263250 16 16            TRUE
## 2  0.2026230 0.1256415 16 16            TRUE
## 3  0.1137269 0.1252021 16 16            TRUE
## 4 -0.1544626 0.1253728 16 16            TRUE
## 5 -0.3400369 0.1449219 14 14           FALSE
## 6  0.3519809 0.1506459 14 13           FALSE
cat("Total number of effect sizes calculated:", nrow(effect_sizes), "\n")
## Total number of effect sizes calculated: 677

5. Overall Meta-Analysis

# Perform a random-effects meta-analysis across all outcomes
overall_model <- rma(yi, vi, data = effect_sizes, method = "REML")
summary(overall_model)
## 
## Random-Effects Model (k = 677; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc   
## -651.6010  1303.2020  1307.2020  1316.2344  1307.2198   
## 
## tau^2 (estimated amount of total heterogeneity): 0.2366 (SE = 0.0189)
## tau (square root of estimated tau^2 value):      0.4864
## I^2 (total heterogeneity / total variability):   71.88%
## H^2 (total variability / sampling variability):  3.56
## 
## Test for Heterogeneity:
## Q(df = 676) = 2188.8347, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub      
##   0.2038  0.0229  8.9031  <.0001  0.1589  0.2487  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Create a forest plot for the overall effect with study names
# If there are too many studies, we'll sample a subset to display
if(nrow(effect_sizes) > 30) {
  # Create a representative subset for visualization
  set.seed(42)  # For reproducibility
  forest_data <- effect_sizes[sample(nrow(effect_sizes), 30), ]
  forest_model <- rma(yi, vi, data = forest_data, method = "REML")
  forest(forest_model, 
         slab = forest_data$study,  # Use actual study names
         main = "Overall Effect of Neuromodulation (Sample of Studies)",
         xlab = "Hedges' g")
} else {
  forest(overall_model, 
         slab = effect_sizes$study,  # Use actual study names
         main = "Overall Effect of Neuromodulation",
         xlab = "Hedges' g")
}

# Create a funnel plot with labels for most extreme points
funnel(overall_model, 
       main = "Funnel Plot: Overall Effects",
       xlab = "Hedges' g")

6. Meta-Analysis by Symptom Domain

6.1 Positive Symptoms

# Filter for positive symptoms
positive_data <- effect_sizes %>% filter(domain == "positive_symptoms")

# Check how many studies we have
cat("Number of effect sizes for positive symptoms:", nrow(positive_data), "\n")
## Number of effect sizes for positive symptoms: 119
# Meta-analysis for positive symptoms if we have enough data
if(nrow(positive_data) >= 3) {
  positive_model <- rma(yi, vi, data = positive_data, method = "REML")
  summary(positive_model)
  
  # Forest plot with study names
  forest(positive_model, 
         slab = positive_data$study,  # Use actual study names
         main = "Effect of Neuromodulation on Positive Symptoms",
         xlab = "Hedges' g")
  
  # Funnel plot to check for publication bias with labels
  funnel(positive_model, 
         main = "Funnel Plot: Positive Symptoms")
  
  # Egger's test for funnel plot asymmetry
  regtest(positive_model, model = "lm")
} else {
  cat("Not enough data for meta-analysis of positive symptoms\n")
}

## 
## Regression Test for Funnel Plot Asymmetry
## 
## Model:     weighted regression with multiplicative dispersion
## Predictor: standard error
## 
## Test for Funnel Plot Asymmetry: t = 0.1304, df = 117, p = 0.8965
## Limit Estimate (as sei -> 0):   b = 0.1580 (CI: -0.2698, 0.5858)

6.2 Negative Symptoms

# Filter for negative symptoms
negative_data <- effect_sizes %>% filter(domain == "negative_symptoms")

# Check how many studies we have
cat("Number of effect sizes for negative symptoms:", nrow(negative_data), "\n")
## Number of effect sizes for negative symptoms: 108
# Meta-analysis for negative symptoms if we have enough data
if(nrow(negative_data) >= 3) {
  negative_model <- rma(yi, vi, data = negative_data, method = "REML")
  summary(negative_model)
  
  # Forest plot with study names
  forest(negative_model, 
         slab = negative_data$study,  # Use actual study names
         main = "Effect of Neuromodulation on Negative Symptoms",
         xlab = "Hedges' g")
  
  # Funnel plot with study labels
  funnel(negative_model, 
         main = "Funnel Plot: Negative Symptoms")
  
  # Egger's test
  regtest(negative_model, model = "lm")
} else {
  cat("Not enough data for meta-analysis of negative symptoms\n")
}

## 
## Regression Test for Funnel Plot Asymmetry
## 
## Model:     weighted regression with multiplicative dispersion
## Predictor: standard error
## 
## Test for Funnel Plot Asymmetry: t = -0.4595, df = 106, p = 0.6468
## Limit Estimate (as sei -> 0):   b =  0.5139 (CI: 0.0997, 0.9281)

6.3 Cognitive Domains

# Identify cognitive domains
cognitive_domains <- c("processing_speed", "attention_vigilance", "working_memory", 
                       "verbal_learning", "visual_learning", "reasoning_problem_solving", 
                       "social_cognition", "global_cognition")

# Loop through each cognitive domain
for(domain_name in cognitive_domains) {
  # Filter data for this domain
  domain_data <- effect_sizes %>% filter(domain == domain_name)
  
  # Check if we have enough data
  if(nrow(domain_data) >= 3) {
    cat("\n\n### Analysis for", domain_name, "###\n")
    cat("Number of effect sizes:", nrow(domain_data), "\n")
    
    # Meta-analysis
    domain_model <- rma(yi, vi, data = domain_data, method = "REML")
    print(summary(domain_model))
    
    # Forest plot with actual study names
    forest(domain_model, 
           slab = domain_data$study,  # Use actual study names
           main = paste("Effect of Neuromodulation on", gsub("_", " ", domain_name)),
           xlab = "Hedges' g")
    
    # Funnel plot with study labels for outliers
    funnel(domain_model, 
           main = paste("Funnel Plot:", gsub("_", " ", domain_name)))
  } else {
    cat("\n\nNot enough data for meta-analysis of", domain_name, "\n")
  }
}
## 
## 
## ### Analysis for processing_speed ###
## Number of effect sizes: 41 
## 
## Random-Effects Model (k = 41; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -16.4346   32.8693   36.8693   40.2470   37.1936   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0152 (SE = 0.0218)
## tau (square root of estimated tau^2 value):      0.1234
## I^2 (total heterogeneity / total variability):   14.33%
## H^2 (total variability / sampling variability):  1.17
## 
## Test for Heterogeneity:
## Q(df = 40) = 45.9103, p-val = 0.2405
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub    
##  -0.0386  0.0520  -0.7425  0.4578  -0.1406  0.0633    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for attention_vigilance ###
## Number of effect sizes: 20 
## 
## Random-Effects Model (k = 20; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -9.8724   19.7448   23.7448   25.6337   24.4948   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0296 (SE = 0.0391)
## tau (square root of estimated tau^2 value):      0.1720
## I^2 (total heterogeneity / total variability):   23.81%
## H^2 (total variability / sampling variability):  1.31
## 
## Test for Heterogeneity:
## Q(df = 19) = 26.4044, p-val = 0.1193
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.0222  0.0812  0.2735  0.7844  -0.1370  0.1814    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for working_memory ###
## Number of effect sizes: 67 
## 
## Random-Effects Model (k = 67; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -65.0222  130.0445  134.0445  138.4238  134.2350   
## 
## tau^2 (estimated amount of total heterogeneity): 0.1419 (SE = 0.0427)
## tau (square root of estimated tau^2 value):      0.3767
## I^2 (total heterogeneity / total variability):   61.25%
## H^2 (total variability / sampling variability):  2.58
## 
## Test for Heterogeneity:
## Q(df = 66) = 179.9725, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.0633  0.0615  1.0292  0.3034  -0.0572  0.1838    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for verbal_learning ###
## Number of effect sizes: 38 
## 
## Random-Effects Model (k = 38; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -24.2909   48.5819   52.5819   55.8037   52.9348   
## 
## tau^2 (estimated amount of total heterogeneity): 0.1218 (SE = 0.0480)
## tau (square root of estimated tau^2 value):      0.3490
## I^2 (total heterogeneity / total variability):   61.93%
## H^2 (total variability / sampling variability):  2.63
## 
## Test for Heterogeneity:
## Q(df = 37) = 92.5419, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.1405  0.0749  1.8761  0.0606  -0.0063  0.2873  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for visual_learning ###
## Number of effect sizes: 34 
## 
## Random-Effects Model (k = 34; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -12.3932   24.7864   28.7864   31.7794   29.1864   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0395 (SE = 0.0268)
## tau (square root of estimated tau^2 value):      0.1987
## I^2 (total heterogeneity / total variability):   36.41%
## H^2 (total variability / sampling variability):  1.57
## 
## Test for Heterogeneity:
## Q(df = 33) = 51.5074, p-val = 0.0211
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.1074  0.0580  1.8508  0.0642  -0.0063  0.2211  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for reasoning_problem_solving ###
## Number of effect sizes: 64 
## 
## Random-Effects Model (k = 64; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -35.0004   70.0009   74.0009   78.2871   74.2009   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0453 (SE = 0.0240)
## tau (square root of estimated tau^2 value):      0.2128
## I^2 (total heterogeneity / total variability):   33.59%
## H^2 (total variability / sampling variability):  1.51
## 
## Test for Heterogeneity:
## Q(df = 63) = 101.8002, p-val = 0.0014
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub    
##  -0.0567  0.0474  -1.1955  0.2319  -0.1497  0.0363    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for social_cognition ###
## Number of effect sizes: 19 
## 
## Random-Effects Model (k = 19; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -13.9787   27.9574   31.9574   33.7381   32.7574   
## 
## tau^2 (estimated amount of total heterogeneity): 0.1919 (SE = 0.0903)
## tau (square root of estimated tau^2 value):      0.4381
## I^2 (total heterogeneity / total variability):   73.31%
## H^2 (total variability / sampling variability):  3.75
## 
## Test for Heterogeneity:
## Q(df = 18) = 62.3132, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.1028  0.1203  0.8544  0.3929  -0.1330  0.3387    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## 
## ### Analysis for global_cognition ###
## Number of effect sizes: 25 
## 
## Random-Effects Model (k = 25; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -8.8515   17.7030   21.7030   24.0592   22.2745   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0439 (SE = 0.0365)
## tau (square root of estimated tau^2 value):      0.2095
## I^2 (total heterogeneity / total variability):   34.61%
## H^2 (total variability / sampling variability):  1.53
## 
## Test for Heterogeneity:
## Q(df = 24) = 33.7531, p-val = 0.0892
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub    
##   0.1843  0.0730  2.5239  0.0116  0.0412  0.3275  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.4 Total Psychopathology

# Filter for total psychopathology
total_psycho_data <- effect_sizes %>% filter(domain == "total_psychopathology")

# Check how many studies we have
cat("Number of effect sizes for total psychopathology:", nrow(total_psycho_data), "\n")
## Number of effect sizes for total psychopathology: 73
# Meta-analysis for total psychopathology if we have enough data
if(nrow(total_psycho_data) >= 3) {
  total_model <- rma(yi, vi, data = total_psycho_data, method = "REML")
  summary(total_model)
  
  # Forest plot with study names
  forest(total_model, 
         slab = total_psycho_data$study,  # Use actual study names
         main = "Effect of Neuromodulation on Total Psychopathology",
         xlab = "Hedges' g")

   
  # Funnel plot with study labels
  funnel(total_model, 
         main = "Funnel Plot: Total Psychopathology")
  
  # Egger's test
  regtest(total_model, model = "lm")
} else {
  cat("Not enough data for meta-analysis of total psychopathology\n")
}

## 
## Regression Test for Funnel Plot Asymmetry
## 
## Model:     weighted regression with multiplicative dispersion
## Predictor: standard error
## 
## Test for Funnel Plot Asymmetry: t = -0.0731, df = 71, p = 0.9419
## Limit Estimate (as sei -> 0):   b =  0.4299 (CI: -0.0481, 0.9078)

7. Summary Table of Effects by Domain

# Create a summary table of all meta-analytic results
summary_table <- data.frame(
  Domain = character(),
  k = integer(),
  Effect_Size = numeric(),
  CI_Lower = numeric(),
  CI_Upper = numeric(),
  p_value = numeric(),
  I_squared = numeric(),
  stringsAsFactors = FALSE
)

# Add overall results
summary_table <- rbind(summary_table, data.frame(
  Domain = "Overall",
  k = overall_model$k,
  Effect_Size = overall_model$b,
  CI_Lower = overall_model$ci.lb,
  CI_Upper = overall_model$ci.ub,
  p_value = overall_model$pval,
  I_squared = (overall_model$tau2 / (overall_model$tau2 + overall_model$vt)) * 100,
  stringsAsFactors = FALSE
))

# Add domain-specific results if available
domains_to_check <- c("positive_symptoms", "negative_symptoms", "total_psychopathology", 
                      cognitive_domains)

for(domain_name in domains_to_check) {
  domain_data <- effect_sizes %>% filter(domain == domain_name)
  
  if(nrow(domain_data) >= 3) {
    domain_model <- rma(yi, vi, data = domain_data, method = "REML")
    
    summary_table <- rbind(summary_table, data.frame(
      Domain = domain_name,
      k = domain_model$k,
      Effect_Size = domain_model$b,
      CI_Lower = domain_model$ci.lb,
      CI_Upper = domain_model$ci.ub,
      p_value = domain_model$pval,
      I_squared = (domain_model$tau2 / (domain_model$tau2 + domain_model$vt)) * 100,
      stringsAsFactors = FALSE
    ))
  }
}

# Format the table nicely
summary_table$Domain <- gsub("_", " ", summary_table$Domain)
summary_table$Domain <- gsub("\\b(\\w)", "\\U\\1", summary_table$Domain, perl = TRUE)

summary_table$Effect_Size <- round(summary_table$Effect_Size, 2)
summary_table$CI_Lower <- round(summary_table$CI_Lower, 2)
summary_table$CI_Upper <- round(summary_table$CI_Upper, 2)
summary_table$p_value <- round(summary_table$p_value, 3)
summary_table$I_squared <- round(summary_table$I_squared, 1)

# Add significance markers
summary_table$Significance <- ifelse(summary_table$p_value < 0.001, "***",
                                   ifelse(summary_table$p_value < 0.01, "**",
                                         ifelse(summary_table$p_value < 0.05, "*", "")))

# Create effect size with CI column
summary_table$Effect_with_CI <- paste0(summary_table$Effect_Size, " [", 
                                     summary_table$CI_Lower, ", ", 
                                     summary_table$CI_Upper, "]", 
                                     summary_table$Significance)

# Simplified table
simplified_table <- summary_table %>%
  select(Domain, k, Effect_with_CI, p_value, I_squared) %>%
  rename(`Number of Studies` = k,
         `Effect Size [95% CI]` = Effect_with_CI,
         `p-value` = p_value,
         `I²` = I_squared)

# Print the table
kable(simplified_table, 
      caption = "Summary of Meta-Analytic Results by Domain") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = FALSE)
Summary of Meta-Analytic Results by Domain
Domain Number of Studies Effect Size [95% CI] p-value I²
intrcpt Overall 677 0.2 [0.16, 0.25]*** 0.000 71.9
intrcpt1 Positive Symptoms 119 0.19 [0.05, 0.33]** 0.007 80.5
intrcpt2 Negative Symptoms 108 0.41 [0.28, 0.54]*** 0.000 74.6
intrcpt3 Total Psychopathology 73 0.41 [0.27, 0.55]*** 0.000 71.8
intrcpt4 Processing Speed 41 -0.04 [-0.14, 0.06] 0.458 14.3
intrcpt5 Attention Vigilance 20 0.02 [-0.14, 0.18] 0.784 23.8
intrcpt6 Working Memory 67 0.06 [-0.06, 0.18] 0.303 61.2
intrcpt7 Verbal Learning 38 0.14 [-0.01, 0.29] 0.061 61.9
intrcpt8 Visual Learning 34 0.11 [-0.01, 0.22] 0.064 36.4
intrcpt9 Reasoning Problem Solving 64 -0.06 [-0.15, 0.04] 0.232 33.6
intrcpt10 Social Cognition 19 0.1 [-0.13, 0.34] 0.393 73.3
intrcpt11 Global Cognition 25 0.18 [0.04, 0.33]* 0.012 34.6

8. Plot of Effect Sizes by Domain

# Create a dataframe for plotting
plot_data <- summary_table %>%
  select(Domain, Effect_Size, CI_Lower, CI_Upper, p_value) %>%
  mutate(Significant = p_value < 0.05,
         Domain = factor(Domain, levels = rev(Domain)))

# Plot
ggplot(plot_data, aes(x = Effect_Size, y = Domain, xmin = CI_Lower, xmax = CI_Upper, 
                      color = Significant)) +
  geom_point(size = 3) +
  geom_errorbarh(height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkgray") +
  scale_color_manual(values = c("FALSE" = "gray50", "TRUE" = "blue")) +
  labs(title = "Effect of Neuromodulation by Symptom Domain",
       x = "Hedges' g Effect Size",
       y = "",
       color = "Statistically\nSignificant") +
  theme_minimal() +
  theme(legend.position = "bottom")