1# WP21 — meta-analytic models in R (metafor + clubSandwich). Input: data/final/effects.csv (one row per harmonised estimate)2# Columns required: es, es_se, family, taxonomy_id, study_id, model_id, country, continent, income_group, data_mid_year,3# property_type, price_type, estimator_class, time_fe, location_fe, pub_type, n_obs, confidence.4suppressPackageStartupMessages({library(metafor); library(clubSandwich)})5args <- commandArgs(trailingOnly = TRUE)6root <- if (length(args) >= 1) args[1] else "."7d <- read.csv(file.path(root, "data/final/effects.csv"), stringsAsFactors = FALSE)8d <- d[is.finite(d$es) & is.finite(d$es_se) & d$es_se > 0, ]9d$vi <- d$es_se^210out_dir <- file.path(root, "results/econometrics"); dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)11set.seed(20261020)1213pooled <- list(); mods <- list(); pub <- list()14for (fam in unique(d$family)) for (tx in unique(d$taxonomy_id[d$family == fam])) {15 s <- d[d$family == fam & d$taxonomy_id == tx, ]16 k_studies <- length(unique(s$study_id)); n_est <- nrow(s)17 if (k_studies < 5 || n_est < 5) next18 # trim extreme estimates (robustness: both versions reported)19 r1 <- tryCatch(rma.uni(yi = es, vi = vi, data = s, method = "REML"), error = function(e) NULL)20 r3 <- tryCatch(rma.mv(yi = es, V = vi, random = ~ 1 | study_id / model_id, data = s, method = "REML"), error = function(e) NULL)21 if (is.null(r1)) next22 rve <- if (!is.null(r3)) tryCatch(coef_test(r3, vcov = "CR2", cluster = s$study_id), error = function(e) NULL) else NULL23 pi1 <- predict(r1)24 row <- data.frame(family = fam, taxonomy_id = tx, k_studies = k_studies, n_estimates = n_est,25 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,26 tau2 = r1$tau2, I2 = r1$I2, H2 = r1$H2, Q = r1$QE, Q_p = r1$QEp,27 est_ml = if (!is.null(r3)) r3$b[1] else NA, se_ml = if (!is.null(r3)) r3$se[1] else NA,28 sigma2_study = if (!is.null(r3)) r3$sigma2[1] else NA, sigma2_model = if (!is.null(r3)) r3$sigma2[2] else NA,29 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,30 median_es = median(s$es), stringsAsFactors = FALSE)31 pooled[[length(pooled) + 1]] <- row32 # publication selection diagnostics (k >= 10 estimates)33 if (n_est >= 10) {34 eg <- tryCatch(regtest(r1, model = "rma", predictor = "sei"), error = function(e) NULL)35 pet <- tryCatch(rma.uni(yi = es, vi = vi, mods = ~ es_se, data = s, method = "REML"), error = function(e) NULL)36 peese <- tryCatch(rma.uni(yi = es, vi = vi, mods = ~ vi, data = s, method = "REML"), error = function(e) NULL)37 sel <- if (n_est >= 30) tryCatch(selmodel(r1, type = "stepfun", steps = c(0.025, 0.05, 0.10)), error = function(e) NULL) else NULL38 pub[[length(pub) + 1]] <- data.frame(family = fam, taxonomy_id = tx, n_estimates = n_est, k_studies = k_studies,39 egger_z = if (!is.null(eg)) eg$zval else NA, egger_p = if (!is.null(eg)) eg$pval else NA,40 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,41 peese_intercept = if (!is.null(peese)) peese$b[1] else NA, peese_se = if (!is.null(peese)) peese$se[1] else NA,42 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)43 }44 # moderator meta-regression (multilevel) when enough studies45 if (k_studies >= 15) {46 s$pub_wp <- as.integer(s$pub_type != "journal"); s$spatial <- as.integer(s$estimator_class %in% c("spatial"))47 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))48 s$rent <- as.integer(s$price_type == "rent"); s$apt <- as.integer(s$property_type %in% c("condo_apartment", "multifamily"))49 s$nonUS <- as.integer(s$country != "United States")50 f <- ~ period + spatial + ml + pub_wp + lnN + rent + apt + nonUS + time_fe + location_fe51 m <- tryCatch(rma.mv(yi = es, V = vi, mods = f, random = ~ 1 | study_id / model_id, data = s, method = "REML"), error = function(e) NULL)52 if (!is.null(m)) {53 ct <- tryCatch(coef_test(m, vcov = "CR2", cluster = s$study_id), error = function(e) NULL)54 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,55 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)56 }57 }58}59write.csv(do.call(rbind, pooled), file.path(out_dir, "pooled_effects.csv"), row.names = FALSE)60if (length(mods)) write.csv(do.call(rbind, mods), file.path(out_dir, "meta_regression.csv"), row.names = FALSE)61if (length(pub)) write.csv(do.call(rbind, pub), file.path(out_dir, "publication_selection.csv"), row.names = FALSE)62cat("pooled families:", length(pooled), " moderator models:", length(mods), " pub-selection:", length(pub), "\n")63