####### R code for: MERF ECE Effects by Decade Analysis
### Created by Emma Hart (erh2169@tc.columbia.edu)
### Last edited July 2025
### This file includes all of the code for the meta-analytic averages presented in the main text and supplement of "Why are Preschool Programs Becoming Less Effective?"


install.packages(c("haven", "metafor", "broom"))
library(haven)
library(metafor)
library(broom)


analysis_5yrs_school<- read_dta("[YOURPATH]/MERF_ece.dta") 
analysis_ece$w<- analysis_ece$se_clust * analysis_ece$se_clust


###################################
#### MAIN MODELS IN MANUSCRIPT ####
###################################

### FIGURE 1 (also reported in Table S2)
# Meta-analytic averages across all outcomes, decade 1 (to 1999)
d1_p <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0))
d1_f1 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0))
d1_f2 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0))
d1_f3 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0))

# Meta-analytic averages across all outcomes, decade 2 (2000 on)
d2_p <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1))
d2_f1 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1))
d2_f2 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1))
d2_f3 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1))


### In-text estimates
# Statistical test for the differences in post-test effect sizes by decade
test_post_diff <- rma.mv(effectsize, w, mods = ~ intyear_cat2000, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1))

# Meta-analytic averages across all outcomes and waves, decade 1 (to 1999)
d1_all <- rma.mv(effectsize, w, random = ~ 1 | studyid/alignedgroup_id,  data= subset(analysis_ece, intyear_cat2000==0))

# Meta-analytic averages across all outcomes and waves, decade 2 (2000 on)
d2_all <- rma.mv(effectsize, w, random = ~ 1 | studyid/alignedgroup_id,  data= subset(analysis_ece, intyear_cat2000==1))

# Statistical test for the differences in all effect sizes (regardless of wave) by decade
test_all_diff <- rma.mv(effectsize, w, mods = ~ intyear_cat2000 + time_dummy2 + time_dummy3 + time_dummy4, random = ~ 1 | studyid/alignedgroup_id,  data= analysis_ece)


#############################
#### SUPPLEMENTAL MODELS ####
#############################

#### TABLE S2 -- META-ANALYTIC ESTIMATES BY MODERATORS:
# Meta-analytic averages across all outcomes and all years
all_p <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1))
all_f1 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2))
all_f2 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3))
all_f3 <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4))

# Meta-analytic averages across cognitive outcomes, ALL years
all_p_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & cognitived==1))
all_f1_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & cognitived==1))
all_f2_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & cognitived==1))
all_f3_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & cognitived==1))

# Meta-analytic averages across cognitive outcomes, decade 1 (to 1999)
d1_p_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0 & cognitived==1))
d1_f1_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0 & cognitived==1))
d1_f2_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0 & cognitived==1))
d1_f3_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0 & cognitived==1))

# Meta-analytic averages across cognitive outcomes, decade 2 (2000 on)
d2_p_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1 & cognitived==1))
d2_f1_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1 & cognitived==1))
d2_f2_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1 & cognitived==1))
d2_f3_c <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1 & cognitived==1))

# Meta-analytic averages across all social-emotional outcomes, ALL years
all_p_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & socd==1))
all_f1_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & socd==1))
all_f2_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & socd==1))
all_f3_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & socd==1))

# Meta-analytic averages across all social-emotional outcomes, decade 1 (to 1999)
d1_p_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0 & socd==1))
d1_f1_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0 & socd==1)) #only 1 study
d1_f2_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0 & socd==1)) #only 1 study
d1_f3_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0 & socd==1)) 

# Meta-analytic averages across all social-emotional outcomes, decade 2 (2000 on)
d2_p_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1 & socd==1))
d2_f1_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1 & socd==1))
d2_f2_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1 & socd==1))
d2_f3_s <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1 & socd==1))


## ONLY NON-CURRICULAR (PROGRAM) INTERVENTIONS:
# Meta-analytic averages across all outcomes, ALL years
all_p_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & curriculard==0))
all_f1_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & curriculard==0))
all_f2_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & curriculard==0))
all_f3_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & curriculard==0))

# Meta-analytic averages across all outcomes, decade 1 (to 1999)
d1_p_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0 & curriculard==0))
d1_f1_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0 & curriculard==0))
d1_f2_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0 & curriculard==0))
d1_f3_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0 & curriculard==0))

# Meta-analytic averages across all outcomes, decade 2 (2000 on)
d2_p_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1 & curriculard==0))
d2_f1_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1 & curriculard==0))
d2_f2_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1 & curriculard==0))
d2_f3_prog <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1 & curriculard==0)) ## only 1 study


## ONLY CURRICULAR INTERVENTIONS:
# Meta-analytic averages across all outcomes, ALL years
all_p_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & curriculard==1))
all_f1_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & curriculard==1))
all_f2_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & curriculard==1))
all_f3_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & curriculard==1))

# Meta-analytic averages across all outcomes, decade 1 (to 1999)
d1_p_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0 & curriculard==1))
d1_f1_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0 & curriculard==1))
# d1_f2_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0 & curriculard==1)) # terminated - too few obs
d1_f3_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0 & curriculard==1))

# Meta-analytic averages across all outcomes, decade 2 (2000 on)
d2_p_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1 & curriculard==1))
d2_f1_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1 & curriculard==1))
d2_f2_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1 & curriculard==1))
d2_f3_cur <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1 & curriculard==1))



## EXPORTING MAIN ESTIMATES AND FIRST SET OF SUPPLEMENTAL ESTIMATES
model_list <- list(d1_p, d1_f1, d1_f2, d1_f3, d2_p, d2_f1, d2_f2, d2_f3, d1_all, d2_all, all_p, all_f1, all_f2, all_f3, all_p_c, all_f1_c, all_f2_c, all_f3_c, d1_p_c, d1_f1_c, d1_f2_c, d1_f3_c,
                   d2_p_c, d2_f1_c, d2_f2_c, d2_f3_c, all_p_s, all_f1_s, all_f2_s, all_f3_s, d1_p_s, d1_f1_s, d1_f2_s, d1_f3_s, d2_p_s, d2_f1_s, d2_f2_s, d2_f3_s, all_p_prog, all_f1_prog, all_f2_prog,
                   all_f3_prog, d1_p_prog, d1_f1_prog, d1_f2_prog, d1_f3_prog, d2_p_prog, d2_f1_prog, d2_f2_prog, d2_f3_prog, all_p_cur, all_f1_cur, all_f2_cur, all_f3_cur, d1_p_cur, d1_f1_cur, d1_f3_cur,
                   d2_p_cur, d2_f1_cur, d2_f2_cur, d2_f3_cur)
model_names <- c("d1_p", "d1_f1", "d1_f2", "d1_f3", "d2_p", "d2_f1", "d2_f2", "d2_f3", "d1_all", "d2_all", "all_p", "all_f1", "all_f2", "all_f3", "all_p_c", "all_f1_c", "all_f2_c", "all_f3_c",
                 "d1_p_c", "d1_f1_c", "d1_f2_c", "d1_f3_c", "d2_p_c", "d2_f1_c", "d2_f2_c", "d2_f3_c", "all_p_s", "all_f1_s", "all_f2_s", "all_f3_s", "d1_p_s", "d1_f1_s", "d1_f2_s", "d1_f3_s",
                 "d2_p_s", "d2_f1_s", "d2_f2_s", "d2_f3_s", "all_p_prog", "all_f1_prog", "all_f2_prog", "all_f3_prog", "d1_p_prog", "d1_f1_prog", "d1_f2_prog", "d1_f3_prog", "d2_p_prog", "d2_f1_prog",
                 "d2_f2_prog", "d2_f3_prog", "all_p_cur", "all_f1_cur", "all_f2_cur", "all_f3_cur", "d1_p_cur", "d1_f1_cur", "d1_f3_cur", "d2_p_cur", "d2_f1_cur", "d2_f2_cur", "d2_f3_cur")

# Create an empty list to store the tidy outputs
tidy_output_list <- list()

# Iterate over each model and save the tidy output
for (i in 1:length(model_list)) {
  tidy_output <- tidy(model_list[[i]])
  
  # Extract the number of observations from the model object
  obs <- model_list[[i]]$k
  
  # Extract the lower CI bound
  ci_lower <- model_list[[i]]$ci.lb
  
  # Extract the upper CI bound
  ci_upper <- model_list[[i]]$ci.ub
  
  # Add the number of observations/ci's as columns in the tidy output
  tidy_output$obs <- obs
  tidy_output$ci_lower <-ci_lower
  tidy_output$ci_upper <-ci_upper
  
  # Add the model name as a column in the tidy output
  tidy_output$model_name <- model_names[i]
  
  tidy_output_list[[i]] <- tidy_output
}

# Combine all tidy outputs into a single data frame
combined_output <- do.call(rbind, tidy_output_list)
write.csv(combined_output, file = "[YOURPATH]/merf_ece1.csv", row.names = FALSE)



#########################################################
#### SUPPLEMENTAL MODELS WITH ROBUST STANDARD ERRORS ####
#########################################################

#### SUPPLEMENTAL TABLE S3 -- ROBUSTNESS CHECK USING ONLY ALIGNED GROUPS
# Meta-analytic averages across all outcomes and all years
all_p_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & aligned_group_with_post==1))
all_f1_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & aligned_group_with_post==1))
all_f2_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & aligned_group_with_post==1))
all_f3_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & aligned_group_with_post==1))

# Meta-analytic averages across all outcomes, decade 1 (to 1999)
d1_p_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==0 & aligned_group_with_post==1))
d1_f1_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==0 & aligned_group_with_post==1))
d1_f2_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==0 & aligned_group_with_post==1))
d1_f3_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==0 & aligned_group_with_post==1))

# Meta-analytic averages across all outcomes, decade 2 (2000 on)
d2_p_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & intyear_cat2000==1 & aligned_group_with_post==1))
d2_f1_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==2 & intyear_cat2000==1 & aligned_group_with_post==1))
d2_f2_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==3 & intyear_cat2000==1 & aligned_group_with_post==1))
d2_f3_a <- rma.mv(effectsize, w, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==4 & intyear_cat2000==1 & aligned_group_with_post==1))

# Statistical test for the differences in post-test effect sizes by decade
test_post_diff_a <- rma.mv(effectsize, w, mods = ~ intyear_cat2000, random = ~ 1 | studyid,  data= subset(analysis_ece, time_cat4==1 & aligned_group_with_post==1))

# Meta-analytic averages across all outcomes and waves, decade 1 (to 1999)
d1_all_a <- rma.mv(effectsize, w, random = ~ 1 | studyid/alignedgroup_id,  data= subset(analysis_ece, intyear_cat2000==0 & aligned_group_with_post==1))

# Meta-analytic averages across all outcomes and waves, decade 2 (2000 on)
d2_all_a <- rma.mv(effectsize, w, random = ~ 1 | studyid/alignedgroup_id,  data= subset(analysis_ece, intyear_cat2000==1 & aligned_group_with_post==1))

# Statistical test for the differences in all effect sizes (regardless of wave) by decade
test_all_diff_a <- rma.mv(effectsize, w, mods = ~ intyear_cat2000 + time_dummy2 + time_dummy3 + time_dummy4, random = ~ 1 | studyid/alignedgroup_id,  data= subset(analysis_ece, aligned_group_with_post==1))


## EXPORTING SECOND SET OF SUPPLEMENTAL ESTIMATES
model_list <- list(all_p_a, all_f1_a, all_f2_a, all_f3_a, d1_p_a, d1_f1_a, d1_f2_a, d1_f3_a, d2_p_a, d2_f1_a, d2_f2_a, d2_f3_a, test_post_diff_a, d1_all_a, d2_all_a, test_all_diff_a)
model_names <- c("all_p_a", "all_f1_a", "all_f2_a", "all_f3_a", "d1_p_a", "d1_f1_a", "d1_f2_a", "d1_f3_a", "d2_p_a", "d2_f1_a", "d2_f2_a", "d2_f3_a", "test_post_diff_a", "d1_all_a", "d2_all_a", "test_all_diff_a")

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) #Robust standard errors (clustered by study)
  
  # 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]/merf_ece2.csv", row.names = FALSE)
