# WP21 — meta-analytic models in R (metafor + clubSandwich). Input: data/final/effects.csv (one row per harmonised estimate)
# Columns required: es, es_se, family, taxonomy_id, study_id, model_id, country, continent, income_group, data_mid_year,
# property_type, price_type, estimator_class, time_fe, location_fe, pub_type, n_obs, confidence.
suppressPackageStartupMessages({library(metafor); library(clubSandwich)})
args <- commandArgs(trailingOnly = TRUE)
root <- if (length(args) >= 1) args[1] else "."
d <- read.csv(file.path(root, "data/final/effects.csv"), stringsAsFactors = FALSE)
d <- d[is.finite(d$es) & is.finite(d$es_se) & d$es_se > 0, ]
d$vi <- d$es_se^2
out_dir <- file.path(root, "results/econometrics"); dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
set.seed(20261020)

pooled <- list(); mods <- list(); pub <- list()
for (fam in unique(d$family)) for (tx in unique(d$taxonomy_id[d$family == fam])) {
  s <- d[d$family == fam & d$taxonomy_id == tx, ]
  k_studies <- length(unique(s$study_id)); n_est <- nrow(s)
  if (k_studies < 5 || n_est < 5) next
  # trim extreme estimates (robustness: both versions reported)
  r1 <- tryCatch(rma.uni(yi = es, vi = vi, data = s, method = "REML"), error = function(e) NULL)
  r3 <- tryCatch(rma.mv(yi = es, V = vi, random = ~ 1 | study_id / model_id, data = s, method = "REML"), error = function(e) NULL)
  if (is.null(r1)) next
  rve <- if (!is.null(r3)) tryCatch(coef_test(r3, vcov = "CR2", cluster = s$study_id), error = function(e) NULL) else NULL
  pi1 <- predict(r1)
  row <- data.frame(family = fam, taxonomy_id = tx, k_studies = k_studies, n_estimates = n_est,
                    est_re = r1$b[1], se_re = r1$se[1], ci_lo = r1$ci.lb, ci_hi = r1$ci.ub, pi_lo = pi1$pi.lb, pi_hi = pi1$pi.ub,
                    tau2 = r1$tau2, I2 = r1$I2, H2 = r1$H2, Q = r1$QE, Q_p = r1$QEp,
                    est_ml = if (!is.null(r3)) r3$b[1] else NA, se_ml = if (!is.null(r3)) r3$se[1] else NA,
                    sigma2_study = if (!is.null(r3)) r3$sigma2[1] else NA, sigma2_model = if (!is.null(r3)) r3$sigma2[2] else NA,
                    se_rve = if (!is.null(rve)) rve$SE[1] else NA, p_rve = if (!is.null(rve)) rve$p_Satt[1] else NA, df_rve = if (!is.null(rve)) rve$df[1] else NA,
                    median_es = median(s$es), stringsAsFactors = FALSE)
  pooled[[length(pooled) + 1]] <- row
  # publication selection diagnostics (k >= 10 estimates)
  if (n_est >= 10) {
    eg <- tryCatch(regtest(r1, model = "rma", predictor = "sei"), error = function(e) NULL)
    pet <- tryCatch(rma.uni(yi = es, vi = vi, mods = ~ es_se, data = s, method = "REML"), error = function(e) NULL)
    peese <- tryCatch(rma.uni(yi = es, vi = vi, mods = ~ vi, data = s, method = "REML"), error = function(e) NULL)
    sel <- if (n_est >= 30) tryCatch(selmodel(r1, type = "stepfun", steps = c(0.025, 0.05, 0.10)), error = function(e) NULL) else NULL
    pub[[length(pub) + 1]] <- data.frame(family = fam, taxonomy_id = tx, n_estimates = n_est, k_studies = k_studies,
      egger_z = if (!is.null(eg)) eg$zval else NA, egger_p = if (!is.null(eg)) eg$pval else NA,
      pet_intercept = if (!is.null(pet)) pet$b[1] else NA, pet_se = if (!is.null(pet)) pet$se[1] else NA, pet_slope = if (!is.null(pet)) pet$b[2] else NA, pet_slope_p = if (!is.null(pet)) pet$pval[2] else NA,
      peese_intercept = if (!is.null(peese)) peese$b[1] else NA, peese_se = if (!is.null(peese)) peese$se[1] else NA,
      selmodel_est = if (!is.null(sel)) sel$b[1] else NA, selmodel_se = if (!is.null(sel)) sel$se[1] else NA, selmodel_lrt_p = if (!is.null(sel)) sel$LRTp else NA, stringsAsFactors = FALSE)
  }
  # moderator meta-regression (multilevel) when enough studies
  if (k_studies >= 15) {
    s$pub_wp <- as.integer(s$pub_type != "journal"); s$spatial <- as.integer(s$estimator_class %in% c("spatial"))
    s$ml <- as.integer(s$estimator_class %in% c("ml")); s$period <- (s$data_mid_year - 2000) / 10; s$lnN <- log(pmax(s$n_obs, 10))
    s$rent <- as.integer(s$price_type == "rent"); s$apt <- as.integer(s$property_type %in% c("condo_apartment", "multifamily"))
    s$nonUS <- as.integer(s$country != "United States")
    f <- ~ period + spatial + ml + pub_wp + lnN + rent + apt + nonUS + time_fe + location_fe
    m <- tryCatch(rma.mv(yi = es, V = vi, mods = f, random = ~ 1 | study_id / model_id, data = s, method = "REML"), error = function(e) NULL)
    if (!is.null(m)) {
      ct <- tryCatch(coef_test(m, vcov = "CR2", cluster = s$study_id), error = function(e) NULL)
      mods[[length(mods) + 1]] <- data.frame(family = fam, taxonomy_id = tx, term = rownames(m$b), est = as.numeric(m$b), se = m$se, p = m$pval,
                                             se_rve = if (!is.null(ct)) ct$SE else NA, p_rve = if (!is.null(ct)) ct$p_Satt else NA, k_studies = k_studies, n_estimates = n_est, stringsAsFactors = FALSE)
    }
  }
}
write.csv(do.call(rbind, pooled), file.path(out_dir, "pooled_effects.csv"), row.names = FALSE)
if (length(mods)) write.csv(do.call(rbind, mods), file.path(out_dir, "meta_regression.csv"), row.names = FALSE)
if (length(pub)) write.csv(do.call(rbind, pub), file.path(out_dir, "publication_selection.csv"), row.names = FALSE)
cat("pooled families:", length(pooled), " moderator models:", length(mods), " pub-selection:", length(pub), "\n")