Preparation of the datasets

In this section, we aggregate the dataframe at the meta-analysis level to obtain the characteristics of participants (e.g., mean age, % male, etc.) and interventions (e.g., mean intensity and frequency) for each meta-analysis included in the review

Dataset 1

This dataset is used in S1 and S2.B

# identify multiple comparisons group as a unique group
df_raw$year <- gsub("b", "", df_raw$year)

df_raw$author_year <- with(df_raw, paste0(df_raw$author, " (", df_raw$year, ")"))
# calculate the sample size of studies with multiple groups (sum) or outcomes (max)
df_N_CCT <- data.table::setDT(df_raw)[, .(x = sum(N_tot)), by =.(First_author_meta, Year_meta, Intervention_general, Intervention_precise, author, year, Outcome_general, multiple_es, id_row)][, .(N_cor = max(x)), by=.(First_author_meta, Year_meta, Intervention_general, Intervention_precise, author, year, Outcome_general)]

# reinsert corrected sample size in the initial dataset
df1 <- dplyr::left_join(df_raw, df_N_CCT, by = c("First_author_meta", "Year_meta", "Intervention_general", "Intervention_precise", "author", "year", "Outcome_general"))

# aggregate dataset at the outcome level
df_agg_outcome <- df1[, lapply(.SD, function(x) weighted.mean(x, N_cor, na.rm = TRUE)), 
              by = c("First_author_meta", "Year_meta", "Intervention_general", "Intervention_precise", "author", "year", "Outcome_general"), 
              .SDcols = c(# participants
                          "MeanAge", "MeanIQ", "PercMale",
                          # intervention - length
                          "LengthWeek", "DurationMonth", 
                          # intervention -implementer
                          "parent_imp", "prof_imp", 
                          # intervention -setting
                          "sett_home", "sett_clin", "sett_class",
                          # RoB
                          "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                          # outcome
                          "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                          "Active_controls")]

df2 <- dplyr::left_join(data.frame(df_agg_outcome), data.frame(df_N_CCT))

# aggregate dataset at the meta-analysis level
df_agg_meta <- data.table::setDT(df2)[, lapply(.SD, function(x)  weighted.mean(x, N_cor, na.rm = TRUE)), 
                by = c("First_author_meta", "Year_meta", "Intervention_general", "Intervention_precise",
                       "Outcome_general"), 
                .SDcols =  c(# participants
                            "MeanAge", "MeanIQ", "PercMale",
                            # intervention - length
                            "LengthWeek", "DurationMonth", 
                            # intervention -implementer
                            "parent_imp", "prof_imp", 
                            # intervention -setting
                            "sett_home", "sett_clin", "sett_class",
                            # RoB
                            "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                            # outcome
                            "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                            "Active_controls")]

# create some new variables to facilitate data presentation
df_agg_meta$factor <- with(df_agg_meta, paste0(First_author_meta, "_", Year_meta, "_", Intervention_precise, "_", Outcome_general))

df_agg_meta$Paper <- with(df_agg_meta, paste0(First_author_meta, " (", Year_meta, ")"))

df_agg_meta$LengthWeek_round <- round(df_agg_meta$LengthWeek)
df_agg_meta$Intensity <- with(df_agg_meta,
      ifelse(LengthWeek < 2, "< 2 h/wk",
        ifelse(LengthWeek < 5, "2-4 h/wk",
          ifelse(LengthWeek < 10, "5-9 h/wk", 
           ifelse(LengthWeek < 20, "10-19 h/wk",
            ifelse(LengthWeek >= 20, ">= 20h/wk", NA_real_))))))

df_agg_meta$Age <- with(df_agg_meta,
      ifelse(MeanAge < 6, "< 6 yo",
        ifelse(MeanAge < 13, "6-12 yo",
          ifelse(MeanAge < 20, "13-19 yo", 
            ifelse(MeanAge >= 20, ">= 20 yo", NA_real_)))))

df_agg_meta$IQ <- with(df_agg_meta,
      ifelse(MeanIQ < 70, "Low (< 70)",
        ifelse(MeanIQ < 80, "Limit (70-79)",
          ifelse(MeanIQ < 120, "Average (80-119)", 
            ifelse(MeanIQ < 130, "High (120-129)", 
              ifelse(MeanIQ >= 130, "Very high (130+)",
                     NA_real_))))))

df_agg_meta$Duration <- with(df_agg_meta,
      ifelse(DurationMonth < 2, "< 2 months",
        ifelse(DurationMonth <= 6, "2-6 months",
        ifelse(DurationMonth <= 12, "7-12 months",
          ifelse(DurationMonth > 12, "> 12 months",
                 NA_real_)))))

df_agg_meta <- dplyr::left_join(df_agg_meta, df_amstar)

# identify mixed age groups in meta-analyses
mixed_df = data.frame(factor = rep(NA, length(unique(df_raw$factor))), min_age = rep(NA, length(unique(df_raw$factor))), max_age = rep(NA, length(unique(df_raw$factor))), diff = rep(NA, length(unique(df_raw$factor))))
row = 1
for (fact in unique(df_raw$factor)) {
  x_i = df_raw[df_raw$factor == fact, ]
  mixed_df$factor[row] = fact
  mixed_df$min_age[row] = min(x_i$MeanAge, na.rm = TRUE)
  mixed_df$max_age[row] = max(x_i$MeanAge, na.rm = TRUE)
  mixed_df$diff[row] = max(x_i$MeanAge, na.rm = TRUE) - min(x_i$MeanAge, na.rm = TRUE)
  row = row + 1
}

df_agg_meta <- dplyr::left_join(df_agg_meta, mixed_df)
df_agg_meta$Age_bounds = paste0("[", round(df_agg_meta$min_age, 1), " - ", round(df_agg_meta$max_age, 1), "]")

This section is dedicated to identify the most recent and rigorous meta-analysis for each combination of intervention and outcome

inter_ok = c("EIBI", "DEV", "NDBI", "SSG", "TECH", "PMI", "TEACCH", "CBT")
outcome_ok = c("Overall ASD symptoms", "Social communication", 
               "Restricted repetitive behaviors", "Adaptive behaviors", "Cognition(IQ)", 
               "Language(Expressive)", "Language(Receptive)", "Disruptive behaviors")

df_agg_meta$IN_meta = 0
df_agg_meta$Year_meta_num = as.numeric(gsub("a|b", "", df_agg_meta$Year_meta))
df_agg_meta$unique_meta = paste0(df_agg_meta$Intervention_precise, "_", df_agg_meta$Outcome_general)
df_agg_meta$row_id_meta = 1:nrow(df_agg_meta)

for (uniq in unique(df_agg_meta$unique_meta)) {
  df_uniq = df_agg_meta[df_agg_meta$unique_meta == uniq ,]
  
  if (unique(df_uniq$Intervention_precise) %in% inter_ok) {
    if(nrow(df_uniq) > 0) {
      for (age in unique(df_uniq$Age)) {
        # subset for a given age group
        df_uniq_age = df_uniq[df_uniq$Age == age, ]
        # select most recent
        if (any(df_uniq_age$Year_meta_num >= 2016) & any(df_uniq_age$Year_meta_num < 2016)) {
          df_uniq_age = df_uniq_age[df_uniq_age$Year_meta >= 2016, ]
        }
        # select highest amstar score
        row_meta = df_uniq_age[which(df_uniq_age$amstar == max(df_uniq_age$amstar)), ]$row_id_meta
        
        if (length(unique(df_agg_meta[df_agg_meta$row_id_meta %in% row_meta, ]$factor)) > 1) {
          df_agg_meta[df_agg_meta$row_id_meta %in% row_meta & 
                        df_agg_meta$Year_meta_num == max(df_agg_meta$Year_meta_num), ]$IN_meta = 1
          print(paste("AMSTAR equality: ", unique(df_agg_meta[df_agg_meta$row_id_meta %in% row_meta, ]$factor)))
        } else {
          df_agg_meta[df_agg_meta$row_id_meta %in% row_meta, ]$IN_meta = 1
        }
      }
    }
  }
}
## [1] "AMSTAR equality:  Reichow_2018_EIBI_Cognition(IQ)"
## [2] "AMSTAR equality:  Rodgers_2021_EIBI_Cognition(IQ)"
## [1] "AMSTAR equality:  Reichow_2018_EIBI_Adaptive behaviors"
## [2] "AMSTAR equality:  Rodgers_2021_EIBI_Adaptive behaviors"
# exclude additional outcomes from main meta analyses
df_agg_meta$IN_meta = ifelse(!(df_agg_meta$Outcome_general %in% outcome_ok), 0, df_agg_meta$IN_meta)

This section is dedicated to identify a unique occurrence of each individual trial

IN_factor = df_agg_meta[which(df_agg_meta$IN_meta == 1), ]$factor
df_raw$Year_meta_num = as.numeric(gsub("a|b", "", df_raw$Year_meta))
df_raw$IN_trial = 0

for (aut in unique(df_raw$author_year)) {
  df_uniq = df_raw[df_raw$author_year == aut, ]
  for (out in unique(df_uniq$Outcome_general)) {
    df_uniq_out = df_uniq[which(df_uniq$Outcome_general == out), ]
      if (length(unique(df_uniq_out$factor)) > 1) {
          row_trial = df_uniq_out[which(df_uniq_out$amstar == max(df_uniq_out$amstar)), ]$row_id
          if (length(unique(with(df_uniq_out[df_uniq_out$row_id %in% row_trial, ],
                                 paste(First_author_meta, Year_meta)))) > 1) {
            df_raw[which(df_raw$row_id %in% row_trial & 
                           df_raw$Year_meta == max(df_uniq_out$Year_meta_num)), ]$IN_trial = 1
          } else {
            df_raw[df_raw$row_id %in% row_trial, ]$IN_trial = 1
          }
      } else {
        df_raw[which(df_raw$row_id %in% df_uniq_out$row_id), ]$IN_trial = 1
      }
  }
}

df_raw$included_trial = df_raw$IN_trial 

Dataset 2

This dataset is used in S3

# subset dataframe to RCTs with low risk of bias regarding blinding of outcome assessors
df_raw_S1 <- subset(df_raw, rob == "low")

# calculate the sample size of studies with multiple groups (sum) or outcomes (max)
df_N_CCT_S1 <- data.table::setDT(df_raw_S1)[, .(x = sum(N_tot)), by =.(First_author_meta, Year_meta, Intervention_general, Intervention_precise, author, year, Outcome_general, multiple_es, id_row)][, .(N_cor = max(x)), by=.(First_author_meta, Year_meta, Intervention_general, Intervention_precise, author, year, Outcome_general)]

# reinsert corrected sample size in the initial dataset
df1_S1 <- dplyr::left_join(df_raw_S1, df_N_CCT_S1, 
                           by = c("First_author_meta", "Year_meta", "Intervention_general",
                                  "Intervention_precise",
                                  "author", "year", "Outcome_general"))

# aggregate dataset at the outcome level
df_agg_outcome_S1 <- df1_S1[, lapply(.SD, function(x) weighted.mean(x, N_cor, na.rm = TRUE)), 
              by = c("First_author_meta", "Year_meta", "Intervention_general", "Intervention_precise",
                     "author", "year", "Outcome_general"), 
              .SDcols = c(# participants
                          "MeanAge", "MeanIQ", "PercMale",
                          # intervention - length
                          "LengthWeek", "DurationMonth", 
                          # intervention -implementer
                          "parent_imp", "prof_imp", 
                          # intervention -setting
                          "sett_home", "sett_clin", "sett_class",
                          # RoB
                          "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                          # outcome
                          "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                          "Active_controls")]

df2_S1 <- dplyr::left_join(data.frame(df_agg_outcome_S1), data.frame(df_N_CCT_S1))

# aggregate dataset at the meta-analysis level
df_agg_meta_S1 <- data.table::setDT(df2_S1)[, lapply(.SD, function(x)  weighted.mean(x, N_cor, na.rm = TRUE)), 
                by = c("First_author_meta", "Year_meta", "Intervention_general", "Intervention_precise",
                       "Outcome_general"), 
                .SDcols =  c(# participants
                            "MeanAge", "MeanIQ", "PercMale",
                            # intervention - length
                            "LengthWeek", "DurationMonth", 
                            # intervention -implementer
                            "parent_imp", "prof_imp", 
                            # intervention -setting
                            "sett_home", "sett_clin", "sett_class",
                            # RoB
                            "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                            # outcome
                            "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                            "Active_controls")]

# create some new variables to facilitate data presentation
df_agg_meta_S1$factor <- with(df_agg_meta_S1, paste0(First_author_meta, "_", Year_meta, "_", Intervention_precise, "_", Outcome_general))

df_agg_meta_S1$Paper <- with(df_agg_meta_S1, paste0(First_author_meta, " (", Year_meta, ")"))

df_agg_meta_S1$LengthWeek_round <- round(df_agg_meta_S1$LengthWeek)
df_agg_meta_S1$Intensity <- with(df_agg_meta_S1,
      ifelse(LengthWeek < 2, "< 2 h/wk",
        ifelse(LengthWeek < 5, "2-4 h/wk",
          ifelse(LengthWeek < 10, "5-9 h/wk", 
           ifelse(LengthWeek < 20, "10-19 h/wk",
            ifelse(LengthWeek >= 20, ">= 20h/wk", NA_real_))))))

df_agg_meta_S1$Age <- with(df_agg_meta_S1,
      ifelse(MeanAge < 6, "< 6 yo",
        ifelse(MeanAge < 13, "6-12 yo",
          ifelse(MeanAge < 20, "13-19 yo", 
            ifelse(MeanAge >= 20, ">= 20 yo", NA_real_)))))

df_agg_meta_S1$IQ <- with(df_agg_meta_S1,
      ifelse(MeanIQ < 70, "Low (< 70)",
        ifelse(MeanIQ < 80, "Limit (70-79)",
          ifelse(MeanIQ < 120, "Average (80-119)", 
            ifelse(MeanIQ < 130, "High (120-129)", 
              ifelse(MeanIQ >= 130, "Very high (130+)",
                     NA_real_))))))

df_agg_meta_S1$Duration <- with(df_agg_meta_S1,
      ifelse(DurationMonth < 2, "< 2 months",
        ifelse(DurationMonth <= 6, "2-6 months",
        ifelse(DurationMonth <= 12, "7-12 months",
          ifelse(DurationMonth > 12, "> 12 months",
                 NA_real_)))))

# identify mixed age groups in meta-analyses
mixed_df_S1 = data.frame(factor = rep(NA, length(unique(df_raw_S1$factor))), min_age = rep(NA, length(unique(df_raw_S1$factor))), max_age = rep(NA, length(unique(df_raw_S1$factor))), diff = rep(NA, length(unique(df_raw_S1$factor))))
row = 1
for (fact in unique(df_raw_S1$factor)) {
  x_i = df_raw_S1[df_raw_S1$factor == fact, ]
  mixed_df_S1$factor[row] = fact
  mixed_df_S1$min_age[row] = min(x_i$MeanAge, na.rm = TRUE)
  mixed_df_S1$max_age[row] = max(x_i$MeanAge, na.rm = TRUE)
  mixed_df_S1$diff[row] = max(x_i$MeanAge, na.rm = TRUE) - min(x_i$MeanAge, na.rm = TRUE)
  row = row + 1
}

df_agg_meta_S1 <- dplyr::left_join(df_agg_meta_S1, mixed_df_S1)
df_agg_meta_S1$Age_bounds = paste0("[", round(df_agg_meta_S1$min_age, 1), " - ", round(df_agg_meta_S1$max_age, 1), "]")

This dataset is used in S2.C

df_raw_S3_t <- subset(df_raw, IN_trial == 1)
Interventions <- stringr::str_split_fixed(df_raw_S3_t$Intervention_all, ' / ', 3)
df_raw_S3_t$intervention1 <- ifelse(Interventions[, 1] == "", 
                                   NA, Interventions[, 1])
df_raw_S3_t$intervention2 <- ifelse(Interventions[, 2] == "", 
                                   NA, Interventions[, 2])
df_raw_S3_t$intervention3 <- ifelse(Interventions[, 3] == "", 
                                   NA, Interventions[, 3])

df_raw_S3 = df_raw_S3_t %>% 
  tidyr::pivot_longer(cols = c(intervention1, intervention2, intervention3),
               names_to = "intervention_rank",
               values_to = "Intervention_long")
df_raw_S3 = subset(df_raw_S3, !is.na(Intervention_long))

df_raw_S3$Age <- with(df_raw_S3,
      ifelse(MeanAge < 6, "< 6 yo",
        ifelse(MeanAge < 13, "6-12 yo",
          ifelse(MeanAge < 20, "13-19 yo", 
            ifelse(MeanAge >= 20, ">= 20 yo", NA_real_)))))

df_raw_S3$factor <- paste0(df_raw_S3$Intervention_long, "_",
                          df_raw_S3$Outcome_general, "_", 
                          df_raw_S3$Age)



# calculate the sample size of studies with multiple groups (sum) or outcomes (max)
df_N_CCT_S3 <- data.table::setDT(df_raw_S3)[, .(x = sum(N_tot)), by =.(Intervention_long, Age, Outcome_general, author, year, multiple_es, id_row)][, .(N_cor = max(x)), by=.(Intervention_long, Age, Outcome_general, author, year)]

# reinsert corrected sample size in the initial dataset
df1_S3 <- dplyr::left_join(df_raw_S3, df_N_CCT_S3, 
                           by = c("Intervention_long",
                                  "Age", "Outcome_general",
                                  "author", "year"))

# aggregate dataset at the outcome level
df_agg_outcome_S3 <- df1_S3[, lapply(.SD, function(x) weighted.mean(x, N_cor, na.rm = TRUE)), 
              by = c("Intervention_long", "Age", "Outcome_general",
                     "author", "year"), 
              .SDcols = c(# participants
                          "MeanAge", "MeanIQ", "PercMale",
                          # intervention - length
                          "LengthWeek", "DurationMonth", 
                          # intervention -implementer
                          "parent_imp", "prof_imp", 
                          # intervention -setting
                          "sett_home", "sett_clin", "sett_class",
                          # RoB
                          "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                          # outcome
                          "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                          "Active_controls")]

df2_S3 <- dplyr::left_join(data.frame(df_agg_outcome_S3), data.frame(df_N_CCT_S3))

# aggregate dataset at the meta-analysis level
df_agg_meta_S3 <- data.table::setDT(df2_S3)[, lapply(.SD, function(x)  weighted.mean(x, N_cor, na.rm = TRUE)), 
              by = c("Intervention_long", "Age", "Outcome_general"), 
                .SDcols =  c(# participants
                            "MeanAge", "MeanIQ", "PercMale",
                            # intervention - length
                            "LengthWeek", "DurationMonth", 
                            # intervention -implementer
                            "parent_imp", "prof_imp", 
                            # intervention -setting
                            "sett_home", "sett_clin", "sett_class",
                            # RoB
                            "Design_num", "RoB_BlindOut_num", "low_RoB_num",
                            # outcome
                            "Outcome_test", "Outcome_questionnaire", "Outcome_observation",
                            "Active_controls")]

# create some new variables to facilitate data presentation
df_agg_meta_S3$factor <- with(df_agg_meta_S3, paste0(Intervention_long, "_", Outcome_general, "_", Age))

df_agg_meta_S3$LengthWeek_round <- round(df_agg_meta_S3$LengthWeek)
df_agg_meta_S3$Intensity <- with(df_agg_meta_S3,
      ifelse(LengthWeek < 2, "< 2 h/wk",
        ifelse(LengthWeek < 5, "2-4 h/wk",
          ifelse(LengthWeek < 10, "5-9 h/wk", 
           ifelse(LengthWeek < 20, "10-19 h/wk",
            ifelse(LengthWeek >= 20, ">= 20h/wk", NA_real_))))))

df_agg_meta_S3$IQ <- with(df_agg_meta_S3,
      ifelse(MeanIQ < 70, "Low (< 70)",
        ifelse(MeanIQ < 80, "Limit (70-79)",
          ifelse(MeanIQ < 120, "Average (80-119)", 
            ifelse(MeanIQ < 130, "High (120-129)", 
              ifelse(MeanIQ >= 130, "Very high (130+)",
                     NA_real_))))))

df_agg_meta_S3$Duration <- with(df_agg_meta_S3,
      ifelse(DurationMonth < 2, "< 2 months",
        ifelse(DurationMonth <= 6, "2-6 months",
        ifelse(DurationMonth <= 12, "7-12 months",
          ifelse(DurationMonth > 12, "> 12 months",
                 NA_real_)))))

# identify mixed age groups in meta-analyses
mixed_df_S3 = data.frame(factor = rep(NA, length(unique(df_raw_S3$factor))), min_age = rep(NA, length(unique(df_raw_S3$factor))), max_age = rep(NA, length(unique(df_raw_S3$factor))), diff = rep(NA, length(unique(df_raw_S3$factor))))
row = 1
for (fact in unique(df_raw_S3$factor)) {
  x_i = df_raw_S3[df_raw_S3$factor == fact, ]
  mixed_df_S3$factor[row] = fact
  mixed_df_S3$min_age[row] = min(x_i$MeanAge, na.rm = TRUE)
  mixed_df_S3$max_age[row] = max(x_i$MeanAge, na.rm = TRUE)
  mixed_df_S3$diff[row] = max(x_i$MeanAge, na.rm = TRUE) - min(x_i$MeanAge, na.rm = TRUE)
  row = row + 1
}

df_agg_meta_S3 <- dplyr::left_join(df_agg_meta_S3, mixed_df_S3)
df_agg_meta_S3$Age_bounds = paste0("[", round(df_agg_meta_S3$min_age, 1), " - ", round(df_agg_meta_S3$max_age, 1), "]")

S1. Main analysis

A. Data analysis

We perform the umbrella review with assessment of the credibility of evidence using algorithmic criteria

df_raw$shared_controls <- NA

umb.main <- metaumbrella::umbrella(df_raw, mult.level = TRUE, method.var = "REML", r = 0.8, verbose = FALSE)
res.transit <- metaumbrella::add.evidence(umb.main, 
    criteria = "Personalized",
    class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                egger_p = .05, pi = "notnull",
                rob = 75, amstar = 4.99),
    class_II = c(total_n = 350, p_value = 1e-3, 
                largest_CI = "notnull", 
                rob = 50, amstar = 3.99),
    class_III = c(total_n = 200, p_value = 1e-3),
    class_IV = c(p_value = 5e-2), verbose = FALSE)

res <- metaumbrella::summary.umbrella(res.transit)
res <- gen_id(res$Factor, res)

We perform some manipulation on the resulting dataset to facilitate data presentation

# merge datasets
res$factor <- with(res, paste0(Author, "_", Year, "_", Intervention, "_", Outcome))
res.main <- dplyr::left_join(res, df_agg_meta)

# estimate the STANDARD ERROR of the pooled ES
ci_lo_transit <- ci_up_transit <- CI_up <- CI_lo <- NA
for (i in 1:nrow(res.main)) {
  ci_lo_transit[i] <- gsub(",.*", "", res.main$value_CI[i])
  ci_up_transit[i] <- gsub(".*, ", "", res.main$value_CI[i])
  CI_lo[i] <- gsub("\\[", "", ci_lo_transit[i])
  CI_up[i] <- gsub("]", "", ci_up_transit[i])
}
CI_lo <- as.numeric(as.character(CI_lo))
CI_up <- as.numeric(as.character(CI_up))
res.main$se <- (CI_up - CI_lo) / 3.92

# Create/recode some variables
res.main$total_n <- round(res.main$total_n)
res.main$Class <- ifelse(res.main$Class == "V", "ns", as.character(res.main$Class))
res.main$Sample_size_studies <- with(res.main, paste0(total_n, " (", n_studies, ")"))
res.main$Effect_size <- paste0("SMD = ", res.main$value, " [", CI_lo, ", ", CI_up, "]")
res.main$Risk_of_bias <- with(res.main, ifelse(rob < 50, "High", 
        ifelse(rob < 75, "Med.", ifelse(rob >= 75, "Low", NA_real_))))
res.main$Active_controls <- paste0(round(res.main$Active_controls), "%")

# Identify main characteristics (setting, outcome type, therapy implementer)
res.main = main_characteristics(res.main)
df_n_out_tot <- do.call(rbind, lapply(split(res.main, res.main$Intervention), n_out))
df_n_out_sub <- do.call(rbind, lapply(split(subset(res.main, IN_meta == 1), subset(res.main, IN_meta == 1)$Intervention), n_out))

df_numb_out_tot <- do.call(rbind, lapply(split(res.main, res.main$Intervention), numb_out_miss))
df_numb_out_sub <- do.call(rbind, lapply(split(subset(res.main, IN_meta == 1), subset(res.main, IN_meta == 1)$Intervention), numb_out_miss))

B. Description of the MA

We provide a general description of the meta-analyses found in the review

############################################
# General description of the meta-analyses #
############################################
paste0(
## Number of eligible articles ---
  paste0("The number of eligible articles identified by our systematic review is equal to: ", length(unique(res.main$Paper)), ". "),
  
  ## Number of meta-analyses and CCTs ---
  paste0("In these articles, a total of ", 
         length(unique(res.main$factor)), 
         " meta-analyses are reported, including ", 
         nrow(unique(cbind(df_raw$author, df_raw$year))) - 1 # lopata/rodgers
         + length(unique(subset(df_raw, multiple_es == "groups")$author)), 
         " unique CCTs and ", nrow(df_raw), " effect size estimates. "),
  
  ## Number of intervention types and outcomes ---
  paste0("The meta-analyses were systematically categorized into 8 intervention types [EIBI (n = ",
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("EIBI", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
         "), NDBI (n = ", 
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("NDBI", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
         "), PMI (n = ", 
         # - 1 for tech pmi
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("PMI", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]) - 1,
          "), SSG (n = ", 
         # - 1 for tech ssg
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("SSG", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]) - 1,
          "), TECH (n = ", 
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("TECH", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
          "), TEACCH (n = ", 
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("TEACCH", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
         "), DEV (n = ", 
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("DEV", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
         "), CBT (n = ", 
         sum(apply(df_n_out_tot, 1, sum)[which(grepl("CBT", names(apply(df_n_out_tot, 1, sum)), fixed = TRUE))]),
         ")] and 8 outcomes [social communication (n = ", 
         sum(df_n_out_tot$SC), 
         "), language (composite: n = ", 
         sum(df_n_out_tot$LT), 
         "; receptive: n = ", 
         sum(df_n_out_tot$LR), 
         "; expressive: n = ", 
         sum(df_n_out_tot$LE), 
          "), IQ (global IQ: n = ", 
         sum(df_n_out_tot$IQ), 
         ", non-verbal IQ: ", 
         sum(df_n_out_tot$nvIQ), 
        "), adaptive behaviors (n = ", 
        sum(df_n_out_tot$AB), 
         "), overall ASD symptoms (n = ", 
        sum(df_n_out_tot$ASD), 
         "), disruptive behaviors (n = ", 
        sum(df_n_out_tot$DIS), 
          "), repetitive and restricted behaviors (n = ", 
        sum(df_n_out_tot$RRB), ")]", ". "),
  
  ## overlapping meta-analyses 
  paste0("After exclusion of overlapping meta-analyses, ", length(unique(subset(res.main, IN_meta == 1)$factor)), " meta-analyses (coming from ", length(unique(subset(res.main, IN_meta == 1)$Paper)), " reports) were included in final analyses. They described ", length(unique(subset(res.main, IN_meta == 1)$factor)), " unique assessements of the efficacy of one of the interventions types on one of the eight possible outcomes.", " ")
)
## [1] "The number of eligible articles identified by our systematic review is equal to: 44. In these articles, a total of 128 meta-analyses are reported, including 192 unique CCTs and 1488 effect size estimates. The meta-analyses were systematically categorized into 8 intervention types [EIBI (n = 36), NDBI (n = 28), PMI (n = 20), SSG (n = 15), TECH (n = 16), TEACCH (n = 7), DEV (n = 4), CBT (n = 2)] and 8 outcomes [social communication (n = 36), language (composite: n = 5; receptive: n = 13; expressive: n = 15), IQ (global IQ: n = 15, non-verbal IQ: 3), adaptive behaviors (n = 15), overall ASD symptoms (n = 16), disruptive behaviors (n = 7), repetitive and restricted behaviors (n = 3)]. After exclusion of overlapping meta-analyses, 46 meta-analyses (coming from 18 reports) were included in final analyses. They described 46 unique assessements of the efficacy of one of the interventions types on one of the eight possible outcomes. "
############################################
##### Description of the meta-analyses #####
############################################
paste0(
## AMSTAR############
paste0("A total of ", 
       nrow(subset(res.main[res.main$IN_meta == 1, ], amstar_rank =="High")), " meta-analyses are ranked 'high' at the AMSTAR 2 tool, ", 
       nrow(subset(res.main[res.main$IN_meta == 1, ], amstar_rank =="Moderate")), " are ranked 'moderate', ",
       nrow(subset(res.main[res.main$IN_meta == 1, ], amstar_rank =="Low")), " are ranked 'low' and ",
       nrow(subset(res.main[res.main$IN_meta == 1, ], amstar_rank =="Critically low")), " are ranked 'critically low'", ". "),

  # Median number of outcomes ---
  paste0("The median number of outcomes assessed for each intervention type was ", median(df_numb_out_sub$N_out), " [min = ", min(df_numb_out_sub$N_out), ", max = ", max(df_numb_out_sub$N_out), "]", ". "),
  
  # Most studied INTERVENTIONS ---
  paste0("The interventions with the highest number of outcomes studied are: ", paste(row.names(subset(df_numb_out_sub, df_numb_out_sub$N_out > 5)), collapse = ', '), ". "),
  paste0("The interventions with the lowest number of outcomes studied are: ", paste(row.names(subset(df_numb_out_sub, df_numb_out_sub$N_out < 4)), collapse = ', '), ". "),
  
  # Least studied OUTCOME ---
  paste0("The two interventions for which RRB was assessed are: ", paste(row.names(df_n_out_sub[df_n_out_sub$RRB==1,]), collapse = ' and '), ". "),
  
  # Median number of studies per meta-analysis  -------------
  paste0("The median number [IQR] of CCTs per meta-analysis is equal to: ", 
         round(median(subset(res.main, IN_meta == 1)$n_studies)), 
         " [", round(median(subset(res.main, IN_meta == 1)$n_studies) - trunc(IQR(subset(res.main, IN_meta == 1)$n_studies))), ", ", 
         round(median(subset(res.main, IN_meta == 1)$n_studies) + IQR(subset(res.main, IN_meta == 1)$n_studies)), "]", ". "),
  
  # N  -------------
  paste0("The median number [IQR] of participants per meta-analysis is equal to: ", 
         round(median(subset(res.main, IN_meta == 1)$total_n)), 
         " [", round(median(subset(res.main, IN_meta == 1)$total_n) - IQR(subset(res.main, IN_meta == 1)$total_n)), ", ",
         round(median(subset(res.main, IN_meta == 1)$total_n) + IQR(subset(res.main, IN_meta == 1)$total_n)), "]", ". "),     
  # RoB  -------------
  paste0("The mean percentage of participants included in studies at low risk of bias per meta-analysis (i.e., included in RCTs with outcome assessors blinded to experimental status) is equal to: ", round(mean(subset(res.main, IN_meta == 1)$rob), 2), "%", ". "),     
  
  # p-value -------------
  paste0("The number of meta-analyses with a nominally significant p-value under the random-effects model is equal to: ", nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02)), " (", round(nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02)) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  paste0("The number of meta-analyses with a p-value < 1e-03 under the random-effects model is equal to: ", nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 1e-03)), " (", round(nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 1e-03)) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),

  paste0("The number of meta-analyses with a p-value < 1e-06 under the random-effects model is equal to: ", nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 1e-06)), " (", round(nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 1e-06)) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  paste0("The number of meta-analyses with a moderate-large SMD (>= .50) is equal to: ", nrow(subset(res.main, IN_meta == 1 & value >= 0.50)), " (", round(nrow(subset(res.main, IN_meta == 1 & value >= 0.50)) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),

  paste0("The number of meta-analyses with a nominally significant p-value OR a moderate-large SMD is equal to: ", nrow(subset(res.main, IN_meta == 1 & (value >= 0.50 | as.numeric(res.main$p_value) < 5e-02))), " (", round(nrow(subset(res.main, IN_meta == 1 & (value >= 0.50 | as.numeric(res.main$p_value) < 5e-02))) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),

  # Largest study -----------
  paste0("Of the ", nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02)), 
         " meta-analyses with statistically significant pooled effect size, ",  
         nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02 & largest_sign == "notnull")), 
         " (", round(nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02 & largest_sign == "notnull")) / nrow(subset(res.main, IN_meta == 1 & as.numeric(res.main$p_value) < 5e-02))*100), 
         "%) were supported by the significant effect of the study with the lowest variance and ", 
  # Prediction interval -----------
nrow(subset(res.main, IN_meta == 1 & PI_sign == "notnull"  & as.numeric(res.main$p_value) < 5e-02)), " (", round(nrow(subset(res.main, IN_meta == 1 & PI_sign == "notnull"  & as.numeric(res.main$p_value) < 5e-02)) / nrow(subset(res.main, IN_meta == 1 & !is.na(PI_sign)  & as.numeric(res.main$p_value) < 5e-02))*100), "%)", " had a prediction interval excluding the null. "),
  
  # inconsistency -----------
  paste0("The number of meta-analyses with a moderate or large inconsistency (I² >= 30%) is equal to: ", nrow(subset(res.main, IN_meta == 1 & I2 >= 30)), " (", round(nrow(subset(res.main, IN_meta == 1 & I2 >= 30)) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  # Small-study effects / publication bias -----------
  paste0("The number of meta-analyses with evidence of small-study effects is equal to: ", nrow(subset(res.main, IN_meta == 1 & (egger_p <= 0.05))), " (", round(nrow(subset(res.main, IN_meta == 1 & (egger_p <= 0.05))) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "))
## [1] "A total of 8 meta-analyses are ranked 'high' at the AMSTAR 2 tool, 0 are ranked 'moderate', 0 are ranked 'low' and 37 are ranked 'critically low'. The median number of outcomes assessed for each intervention type was 5 [min = 1, max = 8]. The interventions with the highest number of outcomes studied are: EIBI, NDBI, PMI. The interventions with the lowest number of outcomes studied are: CBT. The two interventions for which RRB was assessed are: NDBI and SSG. The median number [IQR] of CCTs per meta-analysis is equal to: 6 [2, 10]. The median number [IQR] of participants per meta-analysis is equal to: 238 [68, 409]. The mean percentage of participants included in studies at low risk of bias per meta-analysis (i.e., included in RCTs with outcome assessors blinded to experimental status) is equal to: 19.74%. The number of meta-analyses with a nominally significant p-value under the random-effects model is equal to: 24 (52%). The number of meta-analyses with a p-value < 1e-03 under the random-effects model is equal to: 17 (37%). The number of meta-analyses with a p-value < 1e-06 under the random-effects model is equal to: 6 (13%). The number of meta-analyses with a moderate-large SMD (>= .50) is equal to: 19 (41%). The number of meta-analyses with a nominally significant p-value OR a moderate-large SMD is equal to: 27 (59%). Of the 24 meta-analyses with statistically significant pooled effect size, 13 (54%) were supported by the significant effect of the study with the lowest variance and 9 (41%) had a prediction interval excluding the null. The number of meta-analyses with a moderate or large inconsistency (I² >= 30%) is equal to: 22 (48%). The number of meta-analyses with evidence of small-study effects is equal to: 6 (13%). "
############################################
######## Credibility of evidence### #####
############################################
paste0(
  paste0("The number of meta-analyses associated with a convincing evidence rating is equal to: ", nrow(subset(res.main, IN_meta == 1 & Class == "I")), " (", round(nrow(subset(res.main, IN_meta == 1 & Class == "I")) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  paste0("The number of meta-analyses associated with a highly suggestive evidence rating is equal to: ", nrow(subset(res.main, IN_meta == 1 & Class == "II")), " (", round(nrow(subset(res.main, IN_meta == 1 & Class == "II")) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  
  paste0("The number of meta-analyses associated with a suggestive evidence rating is equal to: ", nrow(subset(res.main, IN_meta == 1 & Class == "III")), " (", round(nrow(subset(res.main, IN_meta == 1 & Class == "III")) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  
  paste0("The number of meta-analyses associated with a weak evidence rating is equal to: ", nrow(subset(res.main, IN_meta == 1 & Class == "IV")), " (", round(nrow(subset(res.main, IN_meta == 1 & Class == "IV")) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". "),
  
  paste0("The number of meta-analyses associated with a non-significant rating is equal to: ", nrow(subset(res.main, IN_meta == 1 & Class == "ns")), " (", round(nrow(subset(res.main, IN_meta == 1 & Class == "ns")) / nrow(subset(res.main, IN_meta == 1))*100), "%)", ". ")
)
## [1] "The number of meta-analyses associated with a convincing evidence rating is equal to: 0 (0%). The number of meta-analyses associated with a highly suggestive evidence rating is equal to: 1 (2%). The number of meta-analyses associated with a suggestive evidence rating is equal to: 10 (22%). The number of meta-analyses associated with a weak evidence rating is equal to: 13 (28%). The number of meta-analyses associated with a non-significant rating is equal to: 22 (48%). "

C. Results

Plot

colnames(res.main)[colnames(res.main) == "Effect_size"] <- "Effect size"
colnames(res.main)[colnames(res.main) == "Sample_size_studies"] <- "Sample size (studies)"
colnames(res.main)[colnames(res.main) == "Type_of_outcome"] <- "Type of outcome"
colnames(res.main)[colnames(res.main) == "Active_controls"] <- "% of active controls"
colnames(res.main)[colnames(res.main) == "Risk_of_bias"] <- "Risk of bias"
res.main$amstar <- paste0(res.main$amstar, "/5")

res.main$Outcome <- rename_outcome(res.main$Outcome)
res.prim.1 <- res.main[, display.cols]
res.prim <- res.prim.1[which(res.prim.1$IN_meta == 1), ]

res.prim$Age <- factor(res.prim$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
res.prim$Class <- factor(res.prim$Class, levels = c("I", "II", "III", "IV", "ns"))
res.prim$Outcome <- factor(res.prim$Outcome,
                                levels = c("Overall ASD sympt.",
                                           "Social-communic.",
                                           "Repetitive behav.",
                                           "Language (Exp)",
                                           "Language (Rec)",
                                           "Adaptive behav.",
                                           "Cognition (IQ)",
                                           "Disruptive behav."))

res.prim <- res.prim[order(res.prim$Outcome,res.prim$Class, res.prim$Age, decreasing = FALSE),]
res.prim <- res.prim[order(res.prim$Age),]

tab.prim <- data.frame(
  Outcome = res.prim$Outcome,
  Intervention = res.prim$Intervention,
  Age = res.prim$Age,
  Class = res.prim$Class,
  RoB = res.prim$'Risk of bias',
  N = round(res.prim$total_n),
  Paper = res.prim$Paper)

forest_modif(x = res.prim[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.prim,
           group = tab.prim$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-2, 2),
           N = res.prim$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.prim[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

Synthesis

Pooled ES

synt.main = subset(res.main, IN_meta == 1)
synt.main$ES_comb = paste(synt.main$value, " ", synt.main$value_CI)
synt.main.ES = synt.main[,c("Intervention", "Outcome", "ES_comb", "Age")] %>% 
                  tidyr::pivot_wider(names_from = Outcome, values_from = ES_comb)


synt.main.ES$Age <- ordered(synt.main.ES$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
synt.main.ES <- synt.main.ES[order(synt.main.ES$Intervention, synt.main.ES$Age), ]
synt.main.ES$Intervention = paste0(synt.main.ES$Intervention, " (", synt.main.ES$Age, ")")
# res.prim$Class <- factor(res.prim$Class, levels = c("I", "II", "III", "IV", "ns"))
# arrange(synt.main.ES, Intervention, Age)

synt.main.ES[is.na(synt.main.ES)] <- "Not assessed"
DT::datatable(subset(synt.main.ES, select = -c(Age)), 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                pageLength = 1600,
                dom = c('t'), 
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '150px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

Selected meta-analyses

synt.main = subset(res.main, IN_meta == 1)
synt.main.ES = synt.main[,c("Intervention", "Outcome", "Paper", "Age")] %>% 
                  tidyr::pivot_wider(names_from = Outcome, values_from = Paper)


synt.main.ES$Age <- ordered(synt.main.ES$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
synt.main.ES <- synt.main.ES[order(synt.main.ES$Intervention, synt.main.ES$Age), ]
synt.main.ES$Intervention = paste0(synt.main.ES$Intervention, " (", synt.main.ES$Age, ")")
# res.prim$Class <- factor(res.prim$Class, levels = c("I", "II", "III", "IV", "ns"))
# arrange(synt.main.ES, Intervention, Age)

synt.main.ES[is.na(synt.main.ES)] <- "Not assessed"
DT::datatable(subset(synt.main.ES, select = -c(Age)), 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                pageLength = 1600,
                dom = c('t'), 
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '150px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

S2. Supplementary analyses

A. Low RoB

This first sensitivity analysis restricts to low risk of bias studies (i.e., RCTs at low risk of detection bias)

Data analysis

df.S1 <- subset(df_raw_S1, factor %in% unique(subset(res.main, IN_meta == 1)$factor))
umb.S1 <- metaumbrella::umbrella(df.S1, mult.level = TRUE, method.var = "REML", r = 0.8, verbose = FALSE)

res.transit.S1 <- metaumbrella::add.evidence(umb.S1, 
    criteria = "Personalised",
    class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                egger_p = .05, pi = "notnull",
                rob = 75, amstar = 4.99),
    class_II = c(total_n = 350, p_value = 1e-3, 
                largest_CI = "notnull", 
                rob = 50, amstar = 3.99),
    class_III = c(total_n = 200, p_value = 1e-3),
    class_IV = c(p_value = 5e-2), verbose = FALSE)

res.S1 <- metaumbrella::summary.umbrella(res.transit.S1)
res.S1 <- gen_id(res.S1$Factor, res.S1)
# merge datasets
res.S1$factor = res.S1$Factor
res.main.S1 <- dplyr::left_join(res.S1, df_agg_meta_S1)

# estimate the STANDARD ERROR of the pooled ES
ci_lo_transit <- ci_up_transit <- CI_up <- CI_lo <- NA
for (i in 1:nrow(res.main.S1)) {
    ci_lo_transit[i] <- gsub(",.*", "", res.main.S1$value_CI[i])
    ci_up_transit[i] <- gsub(".*, ", "", res.main.S1$value_CI[i])
    CI_lo[i] <- gsub("\\[", "", ci_lo_transit[i])
    CI_up[i] <- gsub("]", "", ci_up_transit[i])
}
CI_lo <- as.numeric(as.character(CI_lo))
CI_up <- as.numeric(as.character(CI_up))
res.main.S1$se <- (CI_up - CI_lo) / 3.92

# Create/recode some variables
res.main.S1$total_n <- round(res.main.S1$total_n)
res.main.S1$Class <- ifelse(res.main.S1$Class == "V", "ns", as.character(res.main.S1$Class))
res.main.S1$Sample_size_studies <- with(res.main.S1, paste0(total_n, " (", n_studies, ")"))
res.main.S1$Effect_size <- paste0("SMD = ", res.main.S1$value, " [", CI_lo, ", ", CI_up, "]")
res.main.S1$Risk_of_bias <- with(res.main.S1, ifelse(rob < 50, "High", 
                                               ifelse(rob < 75, "Med.", ifelse(rob >= 75, "Low", NA_real_))))
res.main.S1$Active_controls <- paste0(round(res.main.S1$Active_controls), "%")

# Identify main characteristics (setting, outcome type, therapy implementer)
res.main.S1 = main_characteristics(res.main.S1)

Plot

res.S1.sub <- subset(res.main.S1, n_studies >= 2)
colnames(res.S1.sub)[colnames(res.S1.sub) == "Effect_size"] <- "Effect size"
colnames(res.S1.sub)[colnames(res.S1.sub) == "Sample_size_studies"] <- "Sample size (studies)"
colnames(res.S1.sub)[colnames(res.S1.sub) == "Type_of_outcome"] <- "Type of outcome"
colnames(res.S1.sub)[colnames(res.S1.sub) == "Active_controls"] <- "% of active controls"
colnames(res.S1.sub)[colnames(res.S1.sub) == "Risk_of_bias"] <- "Risk of bias"

res.S1.sub$Outcome <- rename_outcome(res.S1.sub$Outcome)

res.S1.2 <- res.S1.sub[, display.cols[display.cols != "IN_meta"]]

res.S1.2 <- res.S1.2[order(res.S1.2$Outcome,res.S1.2$Class, res.S1.2$Age, decreasing = FALSE),]
res.S1.2 <- res.S1.2[order(res.S1.2$Age),]

res.S1.2$Age <- factor(res.S1.2$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
res.S1.2$Class <- factor(res.S1.2$Class, levels = c("I", "II", "III", "IV", "ns"))
res.S1.2$Outcome <- factor(res.S1.2$Outcome, 
                                levels = c("Overall ASD sympt.", 
                                           "Social-communic.", 
                                           "Repetitive behav.", 
                                           "Language (Exp)", 
                                           "Language (Rec)", 
                                           "Adaptive behav.", 
                                           "Cognition (IQ)", 
                                           "Disruptive behav."))


tab.S1 <- data.frame(
  Outcome = res.S1.2$Outcome,
  Intervention = res.S1.2$Intervention,
  Age = res.S1.2$Age,
  Class = res.S1.2$Class,
  RoB = res.S1.2$'Risk of bias',
  N = round(res.S1.2$total_n),
  Paper = res.S1.2$Paper)

forest_modif(x = res.S1.2[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = F, 
           study_table = tab.S1,
           group = tab.S1$Outcome, 
           type = "study_only",
           text_size = 3.5,
           x_limit = c(-1.1, 1.6),
           N = res.S1.2$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.S1.2, 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))
paste0("The number of meta-analyses with a SMD >= 0.5 is equal to ", nrow(res.prim[res.prim$value >= 0.5,]), " (", round(nrow(res.prim[res.prim$value >= 0.5,]) / nrow(res.prim), 3) * 100, "%)", " in our primary analysis but equal to ", nrow(res.S1.2[res.S1.2$value >= 0.5,]), " (", round(nrow(res.S1.2[res.S1.2$value >= 0.5,]) / nrow(res.S1.2), 3) * 100, "%) in this sensitivity analysis, and the number of meta-analyses with a statistically significant p-value (p<.05) is equal to ", nrow(res.prim[as.numeric(res.prim$p_value) < 0.05,]), " (", round(nrow(res.prim[as.numeric(res.prim$p_value) < 0.05,]) / nrow(res.prim), 3) * 100, "%)", " in our primary analysis but equal to ", nrow(res.S1.2[as.numeric(res.S1.2$p_value) < 0.05,]), " (", round(nrow(res.S1.2[as.numeric(res.S1.2$p_value) < 0.05,]) / nrow(res.S1.2), 3) * 100, "%) in this sensitivity analysis.")
## [1] "The number of meta-analyses with a SMD >= 0.5 is equal to 19 (41.3%) in our primary analysis but equal to 3 (20%) in this sensitivity analysis, and the number of meta-analyses with a statistically significant p-value (p<.05) is equal to 24 (52.2%) in our primary analysis but equal to 4 (26.7%) in this sensitivity analysis."

B. Non-selected MAs

This section presents results of overlapping meta-analyses. To facilitate presentation, each tab contains all meta-analyses for each intervention types

1. EIBI

Plot

res.eibi <- res.main[which(grepl("EIBI", res.main$Intervention, fixed = TRUE)), display.cols]
res.eibi <- res.eibi[order(res.eibi$Class, decreasing = FALSE), ]
res.eibi <- res.eibi[order(res.eibi$Outcome, decreasing = TRUE), ]


tab.eibi <- data.frame(
  Paper = res.eibi$Paper,
  Intervention = res.eibi$Intervention,
  Outcome = res.eibi$Outcome,
  Age = res.eibi$Age,
  Class = res.eibi$Class,
  N = paste0(res.eibi$total_n," "),
  RoB = res.eibi$'Risk of bias',
  Active = res.eibi$`% of active controls`)

forest_modif(x = res.eibi[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.eibi,
           group = tab.eibi$Outcome, 
           type = "study_only",
           text_size = 3.4,
           x_limit = c(-1.7, 5),
           N = res.eibi$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.eibi[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

Discrepancy

df.EIBI.disc <- subset(df_raw, grepl("EIBI", df_raw$Intervention_general, fixed = TRUE) & IN_trial == 1)

df.EIBI.disc$factor <- df.EIBI.disc$Outcome_general
df.EIBI.disc$amstar = min(df.EIBI.disc$amstar)

errors <- metaumbrella::view.errors.umbrella(df.EIBI.disc, return = "data")

df.EIBI.disc$multiple_es[which(df.EIBI.disc$row_id %in% errors[which(errors$column_type_errors == "ERROR"), ]$row_id)] = "outcomes"

umb.EIBI.disc <- metaumbrella::umbrella(df.EIBI.disc, mult.level = TRUE, verbose = FALSE)
## 
res.EIBI.disc = metaumbrella::summary.umbrella(metaumbrella::add.evidence(umb.EIBI.disc, 
                  criteria = "Personalized",
                  class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                              egger_p = .05, pi = "notnull",
                              rob = 75),
                  class_II = c(total_n = 350, p_value = 1e-3, 
                              largest_CI = "notnull", 
                              rob = 50),
                  class_III = c(total_n = 200, p_value = 1e-3),
                  class_IV = c(p_value = 5e-2), verbose = FALSE))

res.EIBI.disc_comp = res.EIBI.disc[res.EIBI.disc$Factor %in% mainsympt, ]
res.EIBI.disc_comp = res.EIBI.disc_comp[order(res.EIBI.disc_comp$Factor), ]
res.EIBI.disc_comp = res.EIBI.disc_comp[, c("Factor", "Class", "measure", "value", "value_CI", "p_value", "n_studies", "total_n")]

res.main.EIBI_comp = res.main[res.main$Outcome_general %in% mainsympt & res.main$Intervention == "EIBI" & res.main$IN_meta == 1, ]
res.main.EIBI_comp = res.main.EIBI_comp[order(res.main.EIBI_comp$Outcome_general), ]
res.main.EIBI_comp = res.main.EIBI_comp[, c("Factor", "Class", "measure", "value", "value_CI", "p_value", "n_studies", "total_n")]
colnames(res.main.EIBI_comp) <- paste0(colnames(res.main.EIBI_comp), "_main")
colnames(res.EIBI.disc_comp) <- paste0(colnames(res.EIBI.disc_comp), "_tot")

DT::datatable(cbind(res.main.EIBI_comp, res.EIBI.disc_comp),
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

2. NDBI

Plot

res.ndbi <- res.main[which(grepl("NDBI", res.main$Intervention, fixed = TRUE)), display.cols]
res.ndbi <- res.ndbi[order(res.ndbi$Class, decreasing = FALSE),]
res.ndbi <- res.ndbi[order(res.ndbi$Outcome, decreasing = TRUE),]

tab.ndbi <- data.frame(
  Paper = res.ndbi$Paper,
  Intervention = res.ndbi$Intervention,
  Outcome = res.ndbi$Outcome,
  Age = res.ndbi$Age,
  Class = res.ndbi$Class,
  N = res.ndbi$total_n,
  RoB = res.ndbi$'Risk of bias')

forest_modif(x = res.ndbi[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.ndbi,
           group = paste(tab.ndbi$Intervention, tab.ndbi$Outcome), 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.ndbi$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.ndbi[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

Discrepancy

df.NDBI.disc <- subset(df_raw, grepl("NDBI", df_raw$Intervention_all, fixed = TRUE) & IN_trial == 1 & MeanAge < 6)

df.NDBI.disc$factor <- df.NDBI.disc$Outcome_general
df.NDBI.disc$amstar = min(df.NDBI.disc$amstar)

errors <- metaumbrella::view.errors.umbrella(df.NDBI.disc, return = "data")

df.NDBI.disc$multiple_es[which(df.NDBI.disc$row_id %in% errors[which(errors$column_type_errors == "ERROR"), ]$row_id)] = "outcomes"

umb.NDBI.disc <- metaumbrella::umbrella(df.NDBI.disc, mult.level = TRUE, verbose = FALSE)
## 
res.NDBI.disc = metaumbrella::summary.umbrella(metaumbrella::add.evidence(umb.NDBI.disc, 
                  criteria = "Personalized",
                  class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                              egger_p = .05, pi = "notnull",
                              rob = 75),
                  class_II = c(total_n = 350, p_value = 1e-3, 
                              largest_CI = "notnull", 
                              rob = 50),
                  class_III = c(total_n = 200, p_value = 1e-3),
                  class_IV = c(p_value = 5e-2), verbose = FALSE))

res.NDBI.disc_comp = res.NDBI.disc[res.NDBI.disc$Factor %in% mainsympt, ]
res.NDBI.disc_comp = res.NDBI.disc_comp[order(res.NDBI.disc_comp$Factor), ]
res.NDBI.disc_comp = res.NDBI.disc_comp[, c("Factor", "Class", "measure", "value", "value_CI", "p_value", "n_studies", "total_n")]

res.main.NDBI_comp = res.main[res.main$Outcome_general %in% mainsympt & res.main$Intervention == "NDBI" & res.main$IN_meta == 1, ]
res.main.NDBI_comp = res.main.NDBI_comp[order(res.main.NDBI_comp$Outcome_general), ]
res.main.NDBI_comp = res.main.NDBI_comp[, c("Factor", "Class", "measure", "value", "value_CI", "p_value", "n_studies", "total_n")]

colnames(res.main.NDBI_comp) <- paste0(colnames(res.main.NDBI_comp), "_main")
colnames(res.NDBI.disc_comp) <- paste0(colnames(res.NDBI.disc_comp), "_tot")

DT::datatable(cbind(res.main.NDBI_comp, res.NDBI.disc_comp), 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

3. DEV

Plot

res.dev <- res.main[which(grepl("DEV", res.main$Intervention)), display.cols]# , ]
res.dev <- res.dev[order(res.dev$Class, decreasing = FALSE),]
res.dev <- res.dev[order(res.dev$Outcome, decreasing = TRUE),]

tab.dev <- data.frame(
  Paper = res.dev$Paper,
  Intervention = res.dev$Intervention,
  Outcome = res.dev$Outcome,
  Age = res.dev$Age,
  Class = res.dev$Class,
  N = res.dev$total_n,
  RoB = res.dev$'Risk of bias')

forest_modif(x = res.dev[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.dev,
           group = tab.dev$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.dev$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.dev[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

4. SSG

Plot

res.ssg <- res.main[which(grepl("SSG", res.main$Intervention)), display.cols]# , ]
res.ssg <- res.ssg[order(res.ssg$Age, decreasing = TRUE),]
res.ssg <- res.ssg[order(res.ssg$Class, decreasing = FALSE),]
res.ssg <- res.ssg[order(res.ssg$Outcome, decreasing = TRUE),]

tab.ssg <- data.frame(
  Paper = res.ssg$Paper,
  Intervention = res.ssg$Intervention,
  Outcome = res.ssg$Outcome,
  Age = res.ssg$Age,
  Class = res.ssg$Class,
  N = res.ssg$total_n,
  RoB = res.ssg$'Risk of bias')

forest_modif(x = res.ssg[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.ssg,
           group = tab.ssg$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.ssg$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.ssg[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

5. PMI

Plot

res.pmi <- res.main[which(grepl("PMI", res.main$Intervention)), display.cols]# , ]
res.pmi <- res.pmi[order(res.pmi$Age, decreasing = FALSE),]
res.pmi <- res.pmi[order(res.pmi$Outcome, decreasing = TRUE),]

tab.pmi <- data.frame(
  Paper = res.pmi$Paper,
  Intervention = res.pmi$Intervention,
  Outcome = res.pmi$Outcome,
  Age = res.pmi$Age,
  Class = res.pmi$Class,
  N = res.pmi$total_n,
  RoB = res.pmi$'Risk of bias')

forest_modif(x = res.pmi[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.pmi,
           group = tab.pmi$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.pmi$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.pmi[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

6. TECH

Plot

res.tech <- res.main[which(grepl("TECH", res.main$Intervention_precise, fixed = TRUE)),
                           display.cols]# , ]
res.tech <- res.tech[order(res.tech$Age, decreasing = TRUE),]
res.tech <- res.tech[order(res.tech$Class, decreasing = FALSE),]
res.tech <- res.tech[order(res.tech$Outcome, decreasing = TRUE),]

tab.tech <- data.frame(
  Paper = res.tech$Paper,
  Intervention = res.tech$Intervention,
  Outcome = res.tech$Outcome,
  Age = paste0(res.tech$Age, " "),
  Class = res.tech$Class,
  N = paste0(round(res.tech$total_n), " "),
  RoB = res.tech$'Risk of bias')

forest_modif(x = res.tech[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.tech,
           group = tab.tech$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.tech$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.tech[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

7. CBT

Plot

res.cbt <- res.main[which(grepl("CBT", res.main$Intervention)), display.cols]# , ]
res.cbt <- res.cbt[order(res.cbt$Class, decreasing = FALSE),]
res.cbt <- res.cbt[order(res.cbt$Outcome, decreasing = TRUE),]

tab.cbt <- data.frame(
  Paper = res.cbt$Paper,
  Intervention = res.cbt$Intervention,
  Outcome = res.cbt$Outcome,
  Age = res.cbt$Age,
  Class = res.cbt$Class,
  N = res.cbt$total_n,
  RoB = res.cbt$'Risk of bias')

forest_modif(x = res.cbt[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.cbt,
           group = tab.cbt$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.cbt$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.cbt[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

8. TEACCH

Plot

res.teacch <- res.main[which(grepl("TEACCH", res.main$Intervention)), display.cols]# , ]
res.teacch <- res.teacch[order(res.teacch$Class, decreasing = FALSE),]
res.teacch <- res.teacch[order(res.teacch$Outcome, decreasing = TRUE),]

tab.teacch <- data.frame(
  Paper = res.teacch$Paper,
  Intervention = res.teacch$Intervention,
  Outcome = res.teacch$Outcome,
  Age = res.teacch$Age,
  Class = res.teacch$Class,
  N = res.teacch$total_n,
  RoB = res.teacch$'Risk of bias')

forest_modif(x = res.teacch[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = FALSE, 
           study_table = tab.teacch,
           group = tab.teacch$Outcome, 
           type = "study_only",
           text_size = 3.6,
           x_limit = c(-1.7, 2.5),
           N = res.teacch$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.teacch[, 1:24], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

C. Precise outcomes

This last supplementary analysis present results of all meta-analyses taken independently (i.e., when a given paper reports several meta-analysis of the same intervention on the same outcome, we now present them individually)

df_all <- df_raw
df_all$factor <- with(df_all, paste0(First_author_meta, "_", Year_meta,  "_", Intervention_precise, "_", Outcome_general, "_", Outcome_precise))

umb.all <- metaumbrella::umbrella(df_all, mult.level = TRUE, method.var = "REML", r = 0.8, verbose = FALSE)
## 
res.transit.all <- metaumbrella::add.evidence(umb.all, 
    criteria = "Personalized",
    class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                egger_p = .05, pi = "notnull",
                rob = 75, amstar = 4.99),
    class_II = c(total_n = 350, p_value = 1e-3, 
                largest_CI = "notnull", 
                rob = 50, amstar = 3.99),
    class_III = c(total_n = 200, p_value = 1e-3),
    class_IV = c(p_value = 5e-2), verbose = FALSE)

res.all <- metaumbrella::summary.umbrella(res.transit.all)
res.all = gen_id_tot(res.all$Factor, res.all)
# merge datasets
res.all$factor <- res.all$Factor
res.all.S2 <- dplyr::left_join(res.all, df_agg_meta)

# estimate the STANDARD ERROR of the pooled ES
ci_lo_transit <- ci_up_transit <- CI_up <- CI_lo <- NA
for (i in 1:nrow(res.all.S2)) {
    ci_lo_transit[i] <- gsub(",.*", "", res.all.S2$value_CI[i])
    ci_up_transit[i] <- gsub(".*, ", "", res.all.S2$value_CI[i])
    CI_lo[i] <- gsub("\\[", "", ci_lo_transit[i])
    CI_up[i] <- gsub("]", "", ci_up_transit[i])
}
CI_lo <- as.numeric(as.character(CI_lo))
CI_up <- as.numeric(as.character(CI_up))
res.all.S2$se <- (CI_up - CI_lo) / 3.92

# Create/recode some variables
res.all.S2$total_n <- round(res.all.S2$total_n)
res.all.S2$Class <- ifelse(res.all.S2$Class == "V", "ns", as.character(res.all.S2$Class))
res.all.S2$Sample_size_studies <- with(res.all.S2, paste0(total_n, " (", n_studies, ")"))
res.all.S2$Effect_size <- paste0("SMD = ", res.all.S2$value, " [", CI_lo, ", ", CI_up, "]")
res.all.S2$Risk_of_bias <- with(res.all.S2, ifelse(rob < 50, "High", 
                                               ifelse(rob < 75, "Med.", ifelse(rob >= 75, "Low", NA_real_))))
res.all.S2$Active_controls <- paste0(round(res.all.S2$Active_controls), "%")

Table

colnames(res.all.S2)[colnames(res.all.S2) == "Effect_size"] <- "Effect size"
colnames(res.all.S2)[colnames(res.all.S2) == "Sample_size_studies"] <- "Sample size (studies)"
colnames(res.all.S2)[colnames(res.all.S2) == "Type_of_outcome"] <- "Type of outcome"
colnames(res.all.S2)[colnames(res.all.S2) == "Active_controls"] <- "% of active controls"
colnames(res.all.S2)[colnames(res.all.S2) == "Risk_of_bias"] <- "Risk of bias"

res.all.S2$Outcome <- rename_outcome(res.all.S2$Outcome)


res.all.S2 <- res.all.S2[order(res.all.S2$Outcome,res.all.S2$Class, res.all.S2$Age, decreasing = FALSE),]
res.all.S2 <- res.all.S2[order(res.all.S2$Age),]

res.all.S2$Age <- factor(res.all.S2$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
res.all.S2$Class <- factor(res.all.S2$Class, levels = c("I", "II", "III", "IV", "ns"))

DT::datatable(res.all.S2[,1:34], 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

S3. Exploratory analysis

This unplanned exploratory analysis aims to limit the inclusion of the same trials in different meta-analysis. More precisely, we assessed the efficacy of each intervention type on each outcome but including all trials identified by all the meta-analyses. To avoid inclusion of the same trial multiple times, we separated the results for each age group using the mean age of the sample in each trial.

Data analysis

for (fact in unique(df_raw_S3$factor)) {
  i = which(df_raw_S3$factor == fact)
  df_raw_S3$amstar[i] = min(df_raw_S3[i, ]$amstar)
}

errors <- metaumbrella::view.errors.umbrella(df_raw_S3, return = "data")

df_raw_S3[df_raw_S3$row_id %in% errors[errors$column_type_errors == "ERROR", ]$row_id , ]$multiple_es = "outcomes"
errors2 <- metaumbrella::view.errors.umbrella(df_raw_S3, return = "data")
umb.tot1 <- metaumbrella::umbrella(df_raw_S3, 
                                   mult.level = TRUE, method.var = "REML", 
                                   r = 0.8, verbose = FALSE)

strat.tot = metaumbrella::summary.umbrella(metaumbrella::add.evidence(umb.tot1, 
                                            criteria = "Personalized",
                                            class_I = c(total_n = 500, p_value = 1e-3, 
                                                        I2 = 50, pi = "notnull",
                                                        egger_p = .05, 
                                                        rob = 75),
                                            class_II = c(total_n = 350, p_value = 1e-3, 
                                                        largest_CI = "notnull", 
                                                        rob = 50),
                                            class_III = c(total_n = 200, p_value = 1e-3),
                                            class_IV = c(p_value = 5e-2), verbose = FALSE))

strat.tot = subset(strat.tot, n_studies > 1)
res.S3 = gen_id_all(strat.tot$Factor, strat.tot)
# merge datasets
res.S3$factor = res.S3$Factor
res.main.S3 <- dplyr::left_join(res.S3, df_agg_meta_S3)

# estimate the STANDARD ERROR of the pooled ES
ci_lo_transit <- ci_up_transit <- CI_up <- CI_lo <- NA
for (i in 1:nrow(res.main.S3)) {
    ci_lo_transit[i] <- gsub(",.*", "", res.main.S3$value_CI[i])
    ci_up_transit[i] <- gsub(".*, ", "", res.main.S3$value_CI[i])
    CI_lo[i] <- gsub("\\[", "", ci_lo_transit[i])
    CI_up[i] <- gsub("]", "", ci_up_transit[i])
}
CI_lo <- as.numeric(as.character(CI_lo))
CI_up <- as.numeric(as.character(CI_up))
res.main.S3$se <- (CI_up - CI_lo) / 3.92

# Create/recode some variables
res.main.S3$total_n <- round(res.main.S3$total_n)
res.main.S3$Class <- ifelse(res.main.S3$Class == "V", "ns", as.character(res.main.S3$Class))
res.main.S3$Sample_size_studies <- with(res.main.S3, paste0(total_n, " (", n_studies, ")"))
res.main.S3$Effect_size <- paste0("SMD = ", res.main.S3$value, " [", CI_lo, ", ", CI_up, "]")
res.main.S3$Risk_of_bias <- with(res.main.S3, ifelse(rob < 50, "High", 
                                               ifelse(rob < 75, "Med.", ifelse(rob >= 75, "Low", NA_real_))))
res.main.S3$Active_controls <- paste0(round(res.main.S3$Active_controls), "%")

# Identify main characteristics (setting, outcome type, therapy implementer)
res.main.S3 = main_characteristics(res.main.S3)

Plot

res.S3.sub <- subset(res.main.S3, n_studies >= 2)
colnames(res.S3.sub)[colnames(res.S3.sub) == "Effect_size"] <- "Effect size"
colnames(res.S3.sub)[colnames(res.S3.sub) == "Sample_size_studies"] <- "Sample size (studies)"
colnames(res.S3.sub)[colnames(res.S3.sub) == "Type_of_outcome"] <- "Type of outcome"
colnames(res.S3.sub)[colnames(res.S3.sub) == "Active_controls"] <- "% of active controls"
colnames(res.S3.sub)[colnames(res.S3.sub) == "Risk_of_bias"] <- "Risk of bias"

res.S3.sub$Outcome <- rename_outcome(res.S3.sub$Outcome)

res.S3.2 <- res.S3.sub[, display.cols.all]

res.S3.2$Age <- factor(res.S3.2$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
res.S3.2$Class <- factor(res.S3.2$Class, levels = c("I", "II", "III", "IV", "ns"))
res.S3.2$Outcome <- factor(res.S3.2$Outcome,
                                levels = c("Overall ASD sympt.",
                                           "Social-communic.",
                                           "Repetitive behav.",
                                           "Language (Exp)",
                                           "Language (Rec)",
                                           "Language",
                                           "Adaptive behav.",
                                           "Disruptive behav.",
                                           "Cognition (IQ)",
                                           "Cognition (n.v. IQ)"))
res.S3.2 <- res.S3.2[order(res.S3.2$Outcome,res.S3.2$Class, res.S3.2$Age, decreasing = FALSE),]
res.S3.2 <- res.S3.2[order(res.S3.2$Age),]



tab.S3 <- data.frame(
  Outcome = res.S3.2$Outcome,
  Intervention = res.S3.2$Intervention,
  Age = res.S3.2$Age,
  Class = res.S3.2$Class,
  RoB = res.S3.2$'Risk of bias',
  N = round(res.S3.2$total_n))

forest_modif(x = res.S3.2[, c("value", "se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = F, 
           study_table = tab.S3,
           group = tab.S3$Outcome, 
           type = "study_only",
           text_size = 3.5,
           x_limit = c(-1.1, 1.6),
           N = res.S3.2$total_n,
           x_breaks = seq(-3, 3, 1))

Table

DT::datatable(res.S3.2, 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                dom = c('ft'), 
                scrollY = "600px", 
                pageLength = 600,
                order = list(0, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '110px',
                       targets = "_all"),
                  list(visible = FALSE,
                       targets = c(22:24)),
                  list(className = 'dt-center', 
                                     targets = "_all"))))

Synthesis

Pooled ES

synt.S3 = res.S3
synt.S3$ES_comb = paste(synt.S3$value, " ", synt.S3$value_CI)
synt.S3.ES = synt.S3[,c("Intervention", "Outcome", "ES_comb", "Age")] %>% 
                  tidyr::pivot_wider(names_from = Outcome, values_from = ES_comb)


synt.S3.ES$Age <- ordered(synt.S3.ES$Age, levels = c("< 6 yo", "6-12 yo", "13-19 yo", ">= 20 yo"))
synt.S3.ES <- synt.S3.ES[order(synt.S3.ES$Intervention, synt.S3.ES$Age), ]
synt.S3.ES$Intervention = paste0(synt.S3.ES$Intervention, " (", synt.S3.ES$Age, ")")
# res.prim$Class <- factor(res.prim$Class, levels = c("I", "II", "III", "IV", "ns"))
# arrange(synt.S3.ES, Intervention, Age)

synt.S3.ES[is.na(synt.S3.ES)] <- "Not assessed"
DT::datatable(synt.S3.ES, 
              rownames = FALSE,
              options = list(  # options
                scrollX = TRUE,
                order = list(1, 'asc'),
                pageLength = 1600,
                dom = c('t'), 
                autoWidth = TRUE,
                columnDefs = list(
                  list(visible = FALSE,
                       targets = 1),
                  list(width = '150px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                                     targets = "_all"))))
##