####### R code for: MERF Soc vs. Cog Manuscript (Main and Supplemental Analyses)
### Created by Emma Hart (erh2169@tc.columbia.edu)
### Last edited 6/28/24
### This file includes all of the code for the meta-regression models in the manuscript and supplement. The regression models in the main manuscript are presented first, followed by tables from the supplement, and descriptive meta-analytic estimates last
### For the regression models, the following naming structure was used to help in lining models up with estimates in the manuscript: tableXXcXXpXX where tableXX indicates the table number ("S" for supplemental models), cXX indicates the column number and pXX indicates the panel
### Note RVE estimation is baked in to the printing syntax (not included model by model); after you save the table with the model output, you'll be able to view the RVE adjustments (denoted with "CI_" in exported table)
### Some models will throw error warnings when the number of post-test effect sizes does not align with the number of follow-up effect sizes being estimated. This is expected.

install.packages("haven")
library(haven)
install.packages("metafor")
library(metafor)
install.packages("sandwich")
library(sandwich)
install.packages("dplyr")
library(dplyr)
install.packages("broom")
library(broom)
install.packages("tidyverse")
library(tidyverse) 
install.packages("kableExtra")
library(kableExtra)
install.packages("robumeta")
library(robumeta) 
install.packages("metafor")
library(metafor)
install.packages("clubSandwich")
library(clubSandwich)
install.packages("weightr")
library(weightr)

#### MAIN REGRESSION MODELS IN MANUSCRIPT

analysis_aligned_viablegroups_wide <- read_dta("[YOURPATH]/MERF_soc_cog_wide.dta") 
dat <- analysis_aligned_viablegroups_wide

### Setting up variance vars 
dat$vi1 <- dat$se1_clust * dat$se1_clust
dat$vi2 <- dat$se2_clust * dat$se2_clust
dat$vi3 <- dat$se3_clust * dat$se3_clust
dat$viz <- dat$sez_clust * dat$sez_clust

dat$w <- 1

### Table 4 (main regression results)
## NULL MODEL HETEREOGENEITY STATS
# regressions without mods (null model, all outcomes)
table4c0p1 <- rma.mv(effectsize2, vi2, random = ~ 1 | studyid, data=dat)
table4c0p1
table4c0p2 <- rma.mv(effectsize3, vi3, random = ~ 1 | studyid, data=dat)
table4c0p2

# regressions without mods (null model, cog outcomes)
table4c00p1 <- rma.mv(effectsize2, vi2, random = ~ 1 | studyid, data= subset(dat, cognitived==1))
table4c00p1
table4c00p2 <- rma.mv(effectsize3, vi3, random = ~ 1 | studyid, data= subset(dat, cognitived==1))
table4c00p2

# regressions without mods (null model, soc outcomes)
table4c000p1 <- rma.mv(effectsize2, vi2, random = ~ 1 | studyid, data= subset(dat, cognitived==0))
table4c000p1
table4c000p2 <- rma.mv(effectsize3, vi3, random = ~ 1 | studyid, data= subset(dat, cognitived==0))
table4c000p2

## REG MODEL ESTIMATES
# regressions with post-test (all outcomes)
table4c1p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
table4c1p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# regressions with post-test and socd (all outcomes)
table4c2p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1 + socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
table4c2p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1 + socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# regressions with post-test, socd, interaction term (all outcomes)
table4c3p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
table4c3p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# regressions with post-test (cog outcomes)
table4c4p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==1))
table4c4p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==1))

# regressions with post-test (soc outcomes)
table4c5p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==0))
table4c5p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==0))


### Table 5 (3 way intervention type interaction regression results)
# regressions with 3-way interaction (cog outcomes)
table5c1p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1 + cog_inputd + soc_inputd + post_cogint + post_socint, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==1))
table5c1p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1 + cog_inputd + soc_inputd + post_cogint + post_socint, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==1))

# regressions with 3-way interaction (soc outcomes)
table5c2p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1 + soc_inputd + post_socint, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==0))
table5c2p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1 + soc_inputd + post_socint, random = ~ effectsize1 | studyid, struct="GEN", data= subset(dat, cognitived==0))



#### SUPPLEMENTAL REGRESSION MODELS

### Table S3 (reg estimates for >2 year follow-up)
## NULL MODEL HETEREOGENEITY STATS
# regressions with post-test, socd, interaction term (all outcomes)
tableS3c0 <- rma.mv(effectsizez, viz,  random = ~ 1 | studyid,  data=dat)
tableS3c0

## REG MODEL ESTIMATES
# regressions with post-test (all outcomes)
tableS3c1 <- rma.mv(effectsizez, viz, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

### Table S4-- Probing Intervention-Related Differences
# econometric fixed effects
dat$fe <- factor(dat$studyid)

tableS4c1p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd + fe, random = NULL, data=dat)
tableS4c1p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd + fe, random = NULL, data=dat)

## age and intensity covariate models
tableS4c2p1<- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd+effectsize1*z_tx_intended_months_i+effectsize1*z_baseline_age_i+tx_intended_months_i_d+baseline_age_i_d, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS4c2p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd+effectsize1*z_tx_intended_months_i+effectsize1*z_baseline_age_i+tx_intended_months_i_d+baseline_age_i_d, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# intensity <200 hrs
tableS4c3p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_200==1))
tableS4c3p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_200==1))

# intensity <100 hrs
tableS4c4p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_100==1))
tableS4c4p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_100==1))

# several intervention characteristic covariate models
tableS4c5p1<- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd+effectsize1*z_ece_i+effectsize1*z_tx_intended_months_split_i+effectsize1*z_tx_parent+effectsize1*z_time_in_school_coarse+effectsize1*z_txyr_2000+effectsize1*z_broad_inputd+effectsize1*z_sample_selection+ece_i_d+tx_intended_months_split_i_d, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS4c5p2<- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd+effectsize1*z_ece_i+effectsize1*z_tx_intended_months_split_i+effectsize1*z_tx_parent+effectsize1*z_time_in_school_coarse+effectsize1*z_txyr_2000+effectsize1*z_broad_inputd+effectsize1*z_sample_selection+ece_i_d+tx_intended_months_split_i_d, random = ~ effectsize1 | studyid, struct="GEN", data=dat)


### Table S5-- Probing Alternate Analytic Approaches and Concerns
## CHE Model

V2 <- vcalc(
  vi = vi2,                                    # sampling variance of effect size estimate
  cluster = studyid,                           # ID for clusters of dependent observations, drawn from the same or partially overlapping samples
  obs = viablegroup,                           # ID for each unique outcome measure within a cluster
  grp1 = studyname_txname,                     # ID for each unique treatment condition within a cluster
  grp2 = cntrl_grp,                            # ID for each unique control condition within a cluster
  w1 = tx_sample_clean2,                       # sample size in treatment group
  w2 = cntrl_sample_clean2,                    # sample size in control group
  data = subset(dat, !is.na(effectsize2)),     # data
  rho = 0.6                                    # assumed correlation between outcomes
)


tableS5c1p1 <- rma.mv(effectsize2 ~ + effectsize1*socd,
                          V = V2,
                          random = list(~ 1+effectsize1|studyid, ~1 |studyname_txname, ~1 |viablegroup),
                          data = subset(dat, !is.na(effectsize2)),
                          sparse = TRUE,
                          struct="GEN") 
tableS5c1p1
tableS5c1p1_rse <- conf_int(tableS5c1p1, vcov = "CR2", p_values = TRUE) # RVE standard errors
tableS5c1p1_rse ## not in output because of different RE structure

V3 <- vcalc(
  vi = vi3,                                   # sampling variance of effect size estimate
  cluster = studyid,                           # ID for clusters of dependent observations, drawn from the same or partially overlapping samples
  obs = viablegroup,                           # ID for each unique outcome measure within a cluster
  grp1 = studyname_txname,                     # ID for each unique treatment condition within a cluster
  grp2 = cntrl_grp,                            # ID for each unique control condition within a cluster
  w1 = tx_sample_clean3,                       # sample size in treatment group
  w2 = cntrl_sample_clean3,                    # sample size in control group
  data = subset(dat, !is.na(effectsize3)),     # data
  rho = 0.6                                    # assumed correlation between outcomes
)


tableS5c1p2 <- rma.mv(effectsize3 ~ + effectsize1*socd,
                      V = V3,
                      random = list(~ 1+effectsize1|studyid, ~1 |studyname_txname, ~1 |viablegroup),
                      data = subset(dat, !is.na(effectsize3)),
                      sparse = TRUE,
                      verbose= TRUE,
                      struct="GEN")
                    
tableS5c1p2
tableS5c1p2_rse <- conf_int(tableS5c1p2, vcov = "CR2", p_values = TRUE) # RVE standard errors
tableS5c1p2_rse ## not in output because of different RE structure

## random effect at meta-analysis level
tableS5c2p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = list(~1+effectsize1| meta_single, ~1+effectsize1 | studyid), struct="GEN", data=dat)
tableS5c2p1_rse <- conf_int(tableS5c2p1, vcov = "CR2", p_values = TRUE) # RVE standard errors
tableS5c2p1_rse ## not in output because of different RE structure
summary(tableS5c2p1)

tableS5c2p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = list(~ 1 +effectsize1| meta_single, ~1+effectsize1 | studyid), struct="GEN", data=dat)
tableS5c2p2_rse <- conf_int(tableS5c2p2, vcov = "CR2", p_values = TRUE) # RVE standard errors
tableS5c2p2_rse ## not in output because of different RE structure
summary(tableS5c2p2)

# cluster covariate
tableS5c3p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd + cluster, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS5c3p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd + cluster, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# no clustered SE adjustment
dat$vi2_noclustadj <- dat$se2 * dat$se2 # set up variance vars that do not use cluster-adjusted SEs
dat$vi3_noclustadj <- dat$se3 * dat$se3

tableS5c4p1 <- rma.mv(effectsize2, vi2_noclustadj, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS5c4p2 <- rma.mv(effectsize3, vi3_noclustadj, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# unweighted
tableS5c5p1 <- rma.mv(effectsize2, vi2, w, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS5c5p2 <- rma.mv(effectsize3, vi3, w, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

# neg post-tests dropped
tableS5c6p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_drop_negposts==1))
tableS5c6p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=subset(dat, analyticsample_drop_negposts==1))

# "estimated" effects dropped
analysis_aligned_viablegroups_dropped_wide <- read_dta("[YOURPATH]/MERF_soc_cog_supp_estdrop.dta")
dat_est <-analysis_aligned_viablegroups_dropped_wide
dat_est$vi2 <- dat_est$se2_clust * dat_est$se2_clust
dat_est$vi3 <- dat_est$se3_clust * dat_est$se3_clust

tableS5c7p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_est)
tableS5c7p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_est)

# models using post-test control group SD (as opposed to concurrent SD)
  # Code to identify the SD ratio (follow-up to post-test SDs) which was then used to adjust follow-up
  # effect sizes for which we did not have post-test SD information to re-calculate effects
  sd_ratio1 <- rma.mv(sd_ratio2, vi2, random = ~ 1 | studyid , data= dat)
  summary(sd_ratio1) ## determine ratio for adjusting effects for which we don't have post-test SD

  sd_ratio1_mods <- rma.mv(sd_ratio2, vi2, mods = ~socd, random = ~ 1 | studyid , data= dat)
  summary(sd_ratio1_mods) ## determine ratio for SOC and COG skills separately

  sd_ratio2 <- rma.mv(sd_ratio3, vi3, random = ~ 1 | studyid , data= dat)
  summary(sd_ratio2) ## determine ratio for adjusting effects for which we don't have post-test SD

  sd_ratio2_mods <- rma.mv(sd_ratio3, vi3,mods = ~socd, random = ~ 1 | studyid , data= dat)
  summary(sd_ratio2_mods) ## determine ratio for SOC and COG skills separately

tableS5c8p1 <- rma.mv(effectsize_sd2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS5c8p2 <- rma.mv(effectsize_sd3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

### Table S7 (reg estimates for PEESE test)
# se estimated on it's own [ONLY REPORTED IN TEXT]
text1_peese <- rma.mv(effectsize1, vi1, mods= ~se1_clust, random = ~ 1 | studyid, data=dat)
text2_peese <- rma.mv(effectsize2, vi2, mods= ~se2_clust, random = ~ 1 | studyid, data=dat)
text3_peese <- rma.mv(effectsize3, vi3, mods= ~se3_clust, random = ~ 1 | studyid, data=dat)
text4_peese <- rma.mv(effectsizez, viz, mods= ~sez_clust, random = ~ 1 | studyid, data=dat)

# se estimated in main regression model with socd and interaction term 
tableS7c1p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd+se2_clust, random = ~ effectsize1 | studyid, struct="GEN", data=dat)
tableS7c1p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd+se3_clust, random = ~ effectsize1 | studyid, struct="GEN", data=dat)

### Table S8 (vevea-hedges)
# Specifying cutpoints and weights to estimate Vevea and Woods (2005):
weightfunct(dat$effectsize1, dat$vi1, steps=c(0.01, 0.05, 0.50, 1.00), weights=c(1, .9, .7, .5))
weightfunct(dat$effectsize2, dat$vi2, steps=c(0.01, 0.05, 0.50, 1.00), weights=c(1, .9, .7, .5))
weightfunct(dat$effectsize3, dat$vi3, steps=c(0.01, 0.05, 0.50, 1.00), weights=c(1, .9, .7, .5))

### Table S9 (reg estimates for alternate grouping approaches)
## regressions with post-test, socd, interaction term with alternate analytic grouping #1 
analysis_aligned_constructgroups_wide <- read_dta("[YOURPATH]/MERF_soc_cog_supp_altgroups1.dta")
dat_A1 <-analysis_aligned_constructgroups_wide
dat_A1$vi2 <- dat_A1$se2_clust * dat_A1$se2_clust
dat_A1$vi3 <- dat_A1$se3_clust * dat_A1$se3_clust

tableS9c1p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_A1)
tableS9c1p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_A1)

## regressions with post-test, socd, interaction term with alternate analytic grouping #2 
viable_cognitivedgroupings_analysis_wide <- read_dta("[YOURPATH]/MERF_soc_cog_supp_altgroups2.dta")
dat_A2 <-viable_cognitivedgroupings_analysis_wide
dat_A2$vi2 <- dat_A2$se2_clust * dat_A2$se2_clust
dat_A2$vi3 <- dat_A2$se3_clust * dat_A2$se3_clust

tableS9c2p1 <- rma.mv(effectsize2, vi2, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_A2)
tableS9c2p2 <- rma.mv(effectsize3, vi3, mods = ~ effectsize1*socd, random = ~ effectsize1 | studyid, struct="GEN", data=dat_A2)


### export tables
model_list <- list(table4c0p1, table4c0p2, table4c00p1, table4c00p2, table4c000p1, table4c000p2, 
                   table4c1p1, table4c1p2, table4c2p1, table4c2p2, table4c3p1, table4c3p2, table4c4p1, 
                   table4c4p2, table4c5p1, table4c5p2, table5c1p1, table5c1p2, table5c2p1, table5c2p2, 
                   tableS3c0, tableS3c1, tableS4c1p1, tableS4c1p2, tableS4c2p1, tableS4c2p2, tableS4c3p1, 
                   tableS4c3p2, tableS4c4p1, tableS4c4p2, tableS4c5p1, tableS4c5p2, tableS5c3p1, tableS5c3p2, 
                   tableS5c4p1, tableS5c4p2, tableS5c5p1, tableS5c5p2, tableS5c6p1, tableS5c6p2, tableS5c7p1,
                   tableS5c7p2, tableS5c8p1, tableS5c8p2, tableS7c1p1, tableS7c1p2, tableS9c1p1, tableS9c1p2, 
                   tableS9c2p1, tableS9c2p2)

model_names <- c("table4c0p1", "table4c0p2", "table4c00p1", "table4c00p2", "table4c000p1", "table4c000p2",
                 "table4c1p1", "table4c1p2", "table4c2p1", "table4c2p2", "table4c3p1", "table4c3p2", "table4c4p1",
                 "table4c4p2", "table4c5p1", "table4c5p2", "table5c1p1", "table5c1p2", "table5c2p1", "table5c2p2",
                 "tableS3c0", "tableS3c1", "tableS4c1p1", "tableS4c1p2", "tableS4c2p1", "tableS4c2p2", "tableS4c3p1",
                 "tableS4c3p2", "tableS4c4p1", "tableS4c4p2", "tableS4c5p1", "tableS4c5p2", "tableS5c3p1", "tableS5c3p2",
                 "tableS5c4p1", "tableS5c4p2", "tableS5c5p1", "tableS5c5p2", "tableS5c6p1", "tableS5c6p2", "tableS5c7p1",
                 "tableS5c7p2", "tableS5c8p1", "tableS5c8p2", "tableS7c1p1", "tableS7c1p2", "tableS9c1p1", "tableS9c1p2",
                 "tableS9c2p1", "tableS9c2p2")

tidy_output_list <- list()

for (i in 1:length(model_list)) {
  
  tidy_output <- tidy(model_list[[i]])
  
  data <- model_list[[i]]$data
  clustering_variable <- data$studyid 
  CI_result <- conf_int(model_list[[i]], vcov = "CR2", cluster = clustering_variable, p_values = TRUE)
  
  # Create data frames for tau2 values
  tau2_values <- data.frame(tau2_1 = model_list[[i]]$tau2[1], tau2_2 = model_list[[i]]$tau2[2])
  
  # Add CI_result to tidy_output
  tidy_output$CI_dv_multilevel <- CI_result
  
  obs <- model_list[[i]]$k
  tidy_output$obs <- obs
  tidy_output$model_name <- model_names[i]
  
  # Join tau2 values with tidy_output
  tidy_output <- cbind(tidy_output, tau2_values)
  
  QE <- model_list[[i]]$QE
  tidy_output$QE <- QE
  QEdf <- model_list[[i]]$QEdf
  tidy_output$QEdf <- QEdf
  
  tidy_output_list[[i]] <- tidy_output
}

### for removing row names so everything exports cohesively
# Convert tibbles to data frames
for (i in 1:length(tidy_output_list)) {
  tidy_output_list[[i]] <- as.data.frame(tidy_output_list[[i]])
}

# Initialize a list to store data frames with unique row names
unique_tidy_output_list <- list()

# Loop through the original list
for (i in 1:length(tidy_output_list)) {
  # Create unique row names for each data frame
  unique_row_names <- make.names(1:nrow(tidy_output_list[[i]]))
  rownames(tidy_output_list[[i]]) <- unique_row_names
  unique_tidy_output_list[[i]] <- tidy_output_list[[i]]
}

# Convert data frames back to tibbles
for (i in 1:length(unique_tidy_output_list)) {
  unique_tidy_output_list[[i]] <- as_tibble(unique_tidy_output_list[[i]])
}

combined_output <- do.call(rbind, tidy_output_list)
write.csv(combined_output, file = "[YOURPATH]/MERFoutput_regression.csv", row.names = FALSE)



### DESCRIPTIVES (TABLE 2) 

dat$w1 <- dat$se1_clust * dat$se1_clust
dat$w2 <- dat$se2_clust * dat$se2_clust
dat$w3 <- dat$se3_clust * dat$se3_clust
dat$w4 <- dat$se4_clust * dat$se4_clust
dat$w5 <- dat$se5_clust * dat$se5_clust
dat$w6 <- dat$se9_clust * dat$se9_clust


### Table 2 weighted estimates
# unweighted avg. post-test ES, outcome types
desc_langu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==9))

desc_mathu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==11))

desc_cogu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==3))

desc_acau <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==13))

desc_achieu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==1))

desc_genu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==15))

desc_subu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==17))

desc_intu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==8))

desc_extu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==6))

desc_crimeu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, construct_cat==4))

desc_intcu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cog_inputd==1))

desc_intsu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, soc_inputd==1))

desc_intbu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, broad_inputd==1))

# weighted avg. post-test ES, outcome types
desc_lang <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==9))

desc_math <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==11))

desc_cog <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==3))

desc_aca <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==13))

desc_achie <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==1))

desc_gen <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==15))

desc_sub <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==17))

desc_int <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==8))

desc_ext <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==6))

desc_crime <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, construct_cat==4))

desc_intc <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cog_inputd==1))

desc_ints <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, soc_inputd==1))

desc_intb <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, broad_inputd==1))


### DESCRIPTIVES (TABLE 3) 

# meta-analytic means with weighting and random effects
desc_a1 <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = dat)
summary(desc_a1)

desc_a2 <- rma.mv(effectsize2, w2, random = ~ 1 | studyid, data = dat)
summary(desc_a2)

desc_a3 <- rma.mv(effectsize3, w3, random = ~ 1 | studyid, data = dat)
summary(desc_a3)

desc_a4 <- rma.mv(effectsize4, w4, random = ~ 1 | studyid, data = dat)
summary(desc_a4)

desc_a5 <- rma.mv(effectsize5, w5, random = ~ 1 | studyid, data = dat)
summary(desc_a5)

desc_a6 <- rma.mv(effectsize9, w6, random = ~ 1 | studyid, data = dat)
summary(desc_a6)

# meta-analytic means with weighting and random effects for cognitived=0 (SOCIAL EMOTIONAL)
desc_s1 <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s1)

desc_s2 <- rma.mv(effectsize2, w2, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s2)

desc_s3 <- rma.mv(effectsize3, w3, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s3)

desc_s4 <- rma.mv(effectsize4, w4, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s4)

desc_s5 <- rma.mv(effectsize5, w5, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s5)


# meta-analytic means with weighting and random effects for cognitived=1 (COGNITIVE)
desc_c1 <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c1)

desc_c2 <- rma.mv(effectsize2, w2, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c2)

desc_c3 <- rma.mv(effectsize3, w3, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c3)

desc_c4 <- rma.mv(effectsize4, w4, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c4)

desc_c5 <- rma.mv(effectsize5, w5, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c5)

desc_c6 <- rma.mv(effectsize9, w6, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c6)


# Post-test meta-analytic means with weighting and random effects
desc_a1p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = dat)
summary(desc_a1p)

desc_a2p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, follow_up2==1))
summary(desc_a2p)

desc_a3p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, follow_up3==1))
#summary(desc_a3p)

desc_a4p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, follow_up4==1))
summary(desc_a4p)

desc_a5p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, follow_up5==1))
#summary(desc_a5p)

desc_a6p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, follow_up6==1))
summary(desc_a6p)


# Post-test meta-analytic means with weighting and random effects for cognitived=0 (SOCIAL EMOTIONAL)
desc_s1p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0))
summary(desc_s1p)

desc_s2p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0 & follow_up2==1))
#summary(desc_s2p)

desc_s3p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0 & follow_up3==1))
#summary(desc_s3p)

desc_s4p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0 & follow_up4==1))
#summary(desc_s4p)

desc_s5p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==0 & follow_up5==1))
#summary(desc_s5p)


# Post-test meta-analytic means with weighting and random effects for cognitived=1 (COGNITIVE)
desc_c1p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1))
summary(desc_c1p)

desc_c2p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1 & follow_up2==1))
summary(desc_c2p)

desc_c3p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1 & follow_up3==1))
summary(desc_c3p)

desc_c4p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1 & follow_up4==1))
summary(desc_c4p)

desc_c5p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1 & follow_up5==1))
summary(desc_c5p)

desc_c6p <- rma.mv(effectsize1, w1, random = ~ 1 | studyid, data = subset(dat, cognitived==1 & follow_up6==1))
summary(desc_c6p)

#Unweighted means
dat$u <- 1

# unweighted: meta-analytic means unweighted
desc_a1u <- rma.mv(effectsize1, w1, w, random = NULL, data = dat)
#summary(desc_a1)

desc_a2u <- rma.mv(effectsize2, w2, w, random = NULL, data = dat)
#summary(desc_a2)

desc_a3u <- rma.mv(effectsize3, w3, w, random = NULL, data = dat)
#summary(desc_a3)

desc_a4u <- rma.mv(effectsize4, w4, w, random = NULL, data = dat)
#summary(desc_a4)

desc_a5u <- rma.mv(effectsize5, w5, w, random = NULL, data = dat)
#summary(desc_a5)

desc_a6u <- rma.mv(effectsize9, w6, w, random = NULL, data = dat)
#summary(desc_a6)

# meta-analytic means unweighted for cognitived=0 (SOCIAL EMOTIONAL)
desc_s1u <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0))
summary(desc_s1u)

desc_s2u <- rma.mv(effectsize2, w2, w, random = NULL, data = subset(dat, cognitived==0))
#summary(desc_s2)

desc_s3u <- rma.mv(effectsize3, w3, w, random = NULL, data = subset(dat, cognitived==0))
#summary(desc_s3)

desc_s4u <- rma.mv(effectsize4, w4, w, random = NULL, data = subset(dat, cognitived==0))
#summary(desc_s4)

desc_s5u <- rma.mv(effectsize5, w5, w, random = NULL, data = subset(dat, cognitived==0))
#summary(desc_s5)

# meta-analytic means unweighted for cognitived=1 (COGNITIVE)
desc_c1u <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1))
summary(desc_c1u)

desc_c2u <- rma.mv(effectsize2, w2, w, random = NULL, data = subset(dat, cognitived==1))
#summary(desc_c2)

desc_c3u <- rma.mv(effectsize3, w3, w, random = NULL, data = subset(dat, cognitived==1))
summary(desc_c3)

desc_c4u <- rma.mv(effectsize4, w4, w, random = NULL, data = subset(dat, cognitived==1))
#summary(desc_c4)

desc_c5u <- rma.mv(effectsize5, w5, w, random = NULL, data = subset(dat, cognitived==1))
#summary(desc_c5)

desc_c6u <- rma.mv(effectsize9, w6, w, random = NULL, data = subset(dat, cognitived==1))
#summary(desc_c6)

# Unweighted Post-test meta-analytic means with weighting and random effects
desc_a1pu <- rma.mv(effectsize1, w1, w, random = NULL,  data = dat)
#summary(desc_a1p)

desc_a2pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, follow_up2==1))
#summary(desc_a2p)

desc_a3pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, follow_up3==1))
#summary(desc_a3p)

desc_a4pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, follow_up4==1))
summary(desc_a4p)

desc_a5pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, follow_up5==1))
#summary(desc_a5p)

desc_a6pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, follow_up6==1))
#summary(desc_a6p)

# Unweighted Post-test meta-analytic means with weighting and random effects for cognitived=0 (SOCIAL EMOTIONAL)
desc_s1pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0))
summary(desc_s1pu)

desc_s2pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0 & follow_up2==1))
#summary(desc_s2p)

desc_s3pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0 & follow_up3==1))
#summary(desc_s3p)

desc_s4pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0 & follow_up4==1))
#summary(desc_s4p)

desc_s5pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==0 & follow_up5==1))
#summary(desc_s5p)

# Unweighted Post-test meta-analytic means with weighting and random effects for cognitived=1 (COGNITIVE)
desc_c1pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1))
summary(desc_c1pu)

desc_c2pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1 & follow_up2==1))
#summary(desc_c2p)

desc_c3pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1 & follow_up3==1))
#summary(desc_c3p)

desc_c4pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1 & follow_up4==1))
#summary(desc_c4p)

desc_c5pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1 & follow_up5==1))
#summary(desc_c5p)

desc_c6pu <- rma.mv(effectsize1, w1, w, random = NULL, data = subset(dat, cognitived==1 & follow_up6==1))
#summary(desc_c6p)


model_list2 <- list(desc_lang, desc_math, desc_cog, desc_aca, desc_achie, desc_gen, desc_sub, desc_int, desc_ext, desc_crime, desc_intc, desc_ints,
                    desc_intb, desc_langu, desc_mathu, desc_cogu, desc_acau, desc_achieu, desc_genu, desc_subu, desc_intu, desc_extu, desc_crimeu, desc_intcu, desc_intsu, desc_intbu,
                    desc_a1, desc_a2, desc_a3, desc_a4, desc_a5, desc_a6, desc_c1, desc_c2, desc_c3, desc_c4, desc_c5, desc_c6, desc_s1, desc_s2, desc_s3, desc_s4, desc_s5,
                    desc_a1p, desc_a2p, desc_a3p, desc_a4p, desc_a5p, desc_a6p, desc_c1p, desc_c2p, desc_c3p, desc_c4p, desc_c5p, desc_c6p, desc_s1p, desc_s2p, desc_s3p, desc_s4p, desc_s5p,
                    desc_a1u, desc_a2u, desc_a3u, desc_a4u, desc_a5u, desc_a6u, desc_c1u, desc_c2u, desc_c3u, desc_c4u, desc_c5u, desc_c6u, desc_s1u, desc_s2u, desc_s3u, desc_s4u, desc_s5u, 
                    desc_a1pu, desc_a2pu, desc_a3pu, desc_a4pu, desc_a5pu, desc_a6pu, desc_c1pu, desc_c2pu, desc_c3pu, desc_c4pu, desc_c5pu, desc_c6pu, desc_s1pu, desc_s2pu, desc_s3pu, desc_s4pu, desc_s5pu)
model_names2 <- c("desc_lang", "desc_math", "desc_cog", "desc_aca", "desc_achie", "desc_gen", "desc_sub", "desc_int", "desc_ext", "desc_crime", "desc_intc", "desc_ints", "desc_intb", "desc_langu", "desc_mathu", 
                  "desc_cogu", "desc_acau", "desc_achieu", "desc_genu", "desc_subu", "desc_intu", "desc_extu", "desc_crimeu", "desc_intcu", "desc_intsu", "desc_intbu", "desc_a1", "desc_a2", "desc_a3", "desc_a4", 
                  "desc_a5", "desc_a6", "desc_c1", "desc_c2", "desc_c3", "desc_c4", "desc_c5", "desc_c6", "desc_s1", "desc_s2", "desc_s3", "desc_s4", "desc_s5", "desc_a1p", "desc_a2p", "desc_a3p", "desc_a4p", 
                  "desc_a5p", "desc_a6p", "desc_c1p", "desc_c2p", "desc_c3p", "desc_c4p", "desc_c5p", "desc_c6p", "desc_s1p", "desc_s2p", "desc_s3p", "desc_s4p", "desc_s5p", "desc_a1u", "desc_a2u", "desc_a3u", 
                  "desc_a4u", "desc_a5u", "desc_a6u", "desc_c1u", "desc_c2u", "desc_c3u", "desc_c4u", "desc_c5u", "desc_c6u", "desc_s1u", "desc_s2u", "desc_s3u", "desc_s4u", "desc_s5u", "desc_a1pu", "desc_a2pu", 
                  "desc_a3pu", "desc_a4pu", "desc_a5pu", "desc_a6pu", "desc_c1pu", "desc_c2pu", "desc_c3pu", "desc_c4pu", "desc_c5pu", "desc_c6pu", "desc_s1pu", "desc_s2pu", "desc_s3pu", "desc_s4pu", "desc_s5pu")

tidy_output_list2 <- list()

for (i in 1:length(model_list2)) {
  
  tidy_output2 <- tidy(model_list2[[i]])
  
  data <- model_list2[[i]]$data
  clustering_variable2 <- data$studyid 
  CI_result2 <- conf_int(model_list2[[i]], vcov = "CR2", cluster = clustering_variable2, p_values = TRUE)
  
  # Add CI_result to tidy_output
  tidy_output2$CI_dv_multilevel <- CI_result2
  
  obs <- model_list2[[i]]$k
  tidy_output2$obs <- obs
  tidy_output2$model_name <- model_names2[i]
  
  tidy_output_list2[[i]] <- tidy_output2}

### for removing row names so everything exports cohesively
# Convert tibbles to data frames
for (i in 1:length(tidy_output_list2)) {
  tidy_output_list2[[i]] <- as.data.frame(tidy_output_list2[[i]])
}

# Initialize a list to store data frames with unique row names
unique_tidy_output_list2 <- list()

# Loop through the original list
for (i in 1:length(tidy_output_list2)) {
  # Create unique row names for each data frame
  unique_row_names2 <- make.names(1:nrow(tidy_output_list2[[i]]))
  rownames(tidy_output_list2[[i]]) <- unique_row_names2
  unique_tidy_output_list2[[i]] <- tidy_output_list2[[i]]
}

# Convert data frames back to tibbles
for (i in 1:length(unique_tidy_output_list2)) {
  unique_tidy_output_list2[[i]] <- as_tibble(unique_tidy_output_list2[[i]])
}

combined_output <- do.call(rbind, tidy_output_list2)
write.csv(combined_output, file = "[YOURPATH]/MERFoutput_descriptives.csv", row.names = FALSE)