diff --git a/benchmarks/R/benchmark_survey_estimators.R b/benchmarks/R/benchmark_survey_estimators.R new file mode 100644 index 00000000..bbd91306 --- /dev/null +++ b/benchmarks/R/benchmark_survey_estimators.R @@ -0,0 +1,649 @@ +#!/usr/bin/env Rscript +# benchmark_survey_estimators.R +# +# Generate golden values for 4 survey-enabled estimators: +# S1: ImputationDiD — control-only WLS regression +# S2: StackedDiD — stacked WLS with Q-weight x survey weight composition +# S3: SunAbraham — interaction-weighted ATT +# S4: TripleDifference — three-way interaction DDD +# +# Usage: Rscript benchmarks/R/benchmark_survey_estimators.R + +suppressPackageStartupMessages({ + library(survey) + library(jsonlite) +}) + +output_dir <- file.path("benchmarks", "data", "synthetic") +if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE) + +cat("=== Survey Estimator Validation Benchmark ===\n\n") + +# =========================================================================== +# DGP 1: Staggered panel (shared by S1, S2, S3) +# =========================================================================== + +generate_staggered_survey_data <- function(seed = 42) { + set.seed(seed) + + # 4 strata with unequal sizes + strata_sizes <- c(25L, 30L, 45L, 50L) # 150 units total + n_strata <- length(strata_sizes) + n_units <- sum(strata_sizes) + + # PSU assignment: 2-3 PSUs per stratum + psu_per_stratum <- c(2L, 2L, 3L, 3L) # 10 PSUs total + fpc_values <- c(100, 150, 200, 250) + + # 5 time periods + T_max <- 5L + + # Build unit-level data + unit_id <- integer(0) + stratum_vec <- integer(0) + psu_vec <- integer(0) + fpc_vec <- numeric(0) + weight_vec <- numeric(0) + first_treat_vec <- numeric(0) + unit_effect_vec <- numeric(0) + + global_psu_id <- 0L + unit_counter <- 0L + + for (h in seq_len(n_strata)) { + n_h <- strata_sizes[h] + n_psu_h <- psu_per_stratum[h] + fpc_h <- fpc_values[h] + w_h <- fpc_h / n_h # inverse selection probability weight + + # Assign units to PSUs within stratum (roughly equal) + psu_assignment <- rep(seq_len(n_psu_h), length.out = n_h) + + for (i in seq_len(n_h)) { + unit_counter <- unit_counter + 1L + + # Staggered treatment: ~30% first_treat=3, ~20% first_treat=4, ~50% never + frac <- i / n_h + if (frac <= 0.3) { + ft <- 3 + } else if (frac <= 0.5) { + ft <- 4 + } else { + ft <- Inf + } + + unit_id <- c(unit_id, unit_counter) + stratum_vec <- c(stratum_vec, h) + psu_vec <- c(psu_vec, global_psu_id + psu_assignment[i]) + fpc_vec <- c(fpc_vec, fpc_h) + weight_vec <- c(weight_vec, w_h) + first_treat_vec <- c(first_treat_vec, ft) + unit_effect_vec <- c(unit_effect_vec, rnorm(1, 0, 2)) + } + global_psu_id <- global_psu_id + n_psu_h + } + + # Expand to panel: 5 periods + n_obs <- n_units * T_max + panel <- data.frame( + unit = rep(unit_id, each = T_max), + period = rep(seq_len(T_max), times = n_units), + stratum = rep(stratum_vec, each = T_max), + psu = rep(psu_vec, each = T_max), + fpc = rep(fpc_vec, each = T_max), + weight = rep(weight_vec, each = T_max), + first_treat = rep(first_treat_vec, each = T_max) + ) + + # Treatment indicator D_it = 1(period >= first_treat & first_treat < Inf) + panel$D_it <- as.integer(panel$period >= panel$first_treat & is.finite(panel$first_treat)) + + # Covariates: TIME-VARYING (unit-level base + period-specific noise) + # Time-invariant covariates are collinear with unit FE, so we need variation + # within units across time for meaningful regression coefficients. + treated_ever <- as.integer(is.finite(rep(first_treat_vec, each = T_max))) + stratum_expanded <- rep(stratum_vec, each = T_max) + panel$x1 <- rnorm(n_obs) + 0.5 * treated_ever + 0.3 * panel$period + panel$x2 <- as.integer(runif(n_obs) < (0.3 + 0.05 * stratum_expanded + 0.05 * panel$period)) + + # Outcome: Y_it = 10 + alpha_i + 5*period + 5*D_it + 2*x1 + 1.5*x2 + noise + ue_expanded <- rep(unit_effect_vec, each = T_max) + noise <- rnorm(n_obs, 0, 1) + panel$outcome <- 10 + ue_expanded + 5 * panel$period + + 5 * panel$D_it + 2 * panel$x1 + 1.5 * panel$x2 + noise + + panel +} + +# =========================================================================== +# DGP 2: DDD cross-section (for S4) +# =========================================================================== + +generate_ddd_survey_data <- function(seed = 42) { + set.seed(seed) + + n <- 200L + n_per_cell <- 25L # 8 cells x 25 = 200 + + # Strata orthogonal to treatment cells + strata_sizes <- c(40L, 45L, 55L, 60L) + psu_per_stratum <- c(2L, 2L, 3L, 3L) + fpc_values <- c(100, 120, 150, 180) + + # Build strata/PSU assignment for all 200 obs (ordered by strata) + stratum_all <- integer(0) + psu_all <- integer(0) + fpc_all <- numeric(0) + weight_all <- numeric(0) + global_psu_id <- 0L + + for (h in seq_along(strata_sizes)) { + n_h <- strata_sizes[h] + n_psu_h <- psu_per_stratum[h] + fpc_h <- fpc_values[h] + w_h <- fpc_h / n_h + + psu_assignment <- rep(seq_len(n_psu_h), length.out = n_h) + stratum_all <- c(stratum_all, rep(h, n_h)) + psu_all <- c(psu_all, global_psu_id + psu_assignment) + fpc_all <- c(fpc_all, rep(fpc_h, n_h)) + weight_all <- c(weight_all, rep(w_h, n_h)) + global_psu_id <- global_psu_id + n_psu_h + } + + # Build treatment cells: G x P x T, 25 per cell + cells <- expand.grid(G = 0:1, P = 0:1, Ti = 0:1) + group_vec <- integer(0) + partition_vec <- integer(0) + time_vec <- integer(0) + for (r in seq_len(nrow(cells))) { + group_vec <- c(group_vec, rep(cells$G[r], n_per_cell)) + partition_vec <- c(partition_vec, rep(cells$P[r], n_per_cell)) + time_vec <- c(time_vec, rep(cells$Ti[r], n_per_cell)) + } + + # Shuffle to break any ordering correlation, then assign strata + shuffle_idx <- sample(n) + group_vec <- group_vec[shuffle_idx] + partition_vec <- partition_vec[shuffle_idx] + time_vec <- time_vec[shuffle_idx] + + # Outcome: Y = 10 + 2*G + 1.5*P + 3*T + 1*G*T + 1*P*T + 0.5*G*P + 3*G*P*T + noise + noise <- rnorm(n, 0, 1) + outcome <- 10 + 2 * group_vec + 1.5 * partition_vec + 3 * time_vec + + 1 * group_vec * time_vec + 1 * partition_vec * time_vec + + 0.5 * group_vec * partition_vec + + 3 * group_vec * partition_vec * time_vec + noise + + data.frame( + obs_id = seq_len(n), + group = group_vec, + partition = partition_vec, + time = time_vec, + outcome = outcome, + stratum = stratum_all, + psu = psu_all, + fpc = fpc_all, + weight = weight_all + ) +} + +# =========================================================================== +# Helper: extract svyglm results for a named coefficient +# =========================================================================== + +extract_svyglm_results <- function(model, design, coef_name) { + coefs <- coef(model) + ses <- SE(model) + idx <- which(names(coefs) == coef_name) + if (length(idx) == 0) stop(paste0("Coefficient '", coef_name, "' not found in model")) + + att <- as.numeric(coefs[idx]) + se <- as.numeric(ses[idx]) + t_stat <- att / se + df <- degf(design) + + t_crit <- qt(0.975, df) + ci_lower <- att - t_crit * se + ci_upper <- att + t_crit * se + + list(att = att, se = se, t_stat = t_stat, df = df, + ci_lower = ci_lower, ci_upper = ci_upper) +} + +# =========================================================================== +# Generate data +# =========================================================================== + +cat("Generating staggered panel data...\n") +staggered_data <- generate_staggered_survey_data(seed = 42) +cat(sprintf(" %d obs, %d units, %d periods\n", + nrow(staggered_data), length(unique(staggered_data$unit)), + length(unique(staggered_data$period)))) +cat(sprintf(" Cohorts: first_treat=3: %d units, first_treat=4: %d units, never: %d units\n", + length(unique(staggered_data$unit[staggered_data$first_treat == 3])), + length(unique(staggered_data$unit[staggered_data$first_treat == 4])), + length(unique(staggered_data$unit[is.infinite(staggered_data$first_treat)])))) + +cat("\nGenerating DDD data...\n") +ddd_data <- generate_ddd_survey_data(seed = 42) +cat(sprintf(" %d obs\n", nrow(ddd_data))) + +results <- list() + +# =========================================================================== +# S1: ImputationDiD — control-only WLS regression +# =========================================================================== + +cat("\n--- S1: ImputationDiD control-only regression ---\n") + +# Omega_0: all untreated observations +# Never-treated at all periods + pre-treatment obs of eventually-treated units +omega_0 <- staggered_data[is.infinite(staggered_data$first_treat) | + staggered_data$period < staggered_data$first_treat, ] +cat(sprintf(" Omega_0: %d observations\n", nrow(omega_0))) + +design_s1 <- svydesign(ids = ~psu, strata = ~stratum, fpc = ~fpc, + weights = ~weight, data = omega_0) + +fit_s1 <- svyglm(outcome ~ factor(unit) + factor(period) + x1 + x2, + design = design_s1) + +# Extract covariate coefficients +coefs_s1 <- coef(fit_s1) +ses_s1 <- SE(fit_s1) +df_s1 <- degf(design_s1) + +results$s1_imputation_did <- list( + scenario = "S1", + description = "ImputationDiD: control-only WLS regression on Omega_0", + estimator = "survey::svyglm", + coef_x1 = as.numeric(coefs_s1["x1"]), + se_x1 = as.numeric(ses_s1["x1"]), + coef_x2 = as.numeric(coefs_s1["x2"]), + se_x2 = as.numeric(ses_s1["x2"]), + df = df_s1, + n_obs_omega0 = nrow(omega_0) +) + +cat(sprintf(" coef_x1 = %.6f (SE = %.6f)\n", + results$s1_imputation_did$coef_x1, results$s1_imputation_did$se_x1)) +cat(sprintf(" coef_x2 = %.6f (SE = %.6f)\n", + results$s1_imputation_did$coef_x2, results$s1_imputation_did$se_x2)) +cat(sprintf(" df = %d\n", df_s1)) + +# =========================================================================== +# S2: StackedDiD — stacked WLS with Q-weight x survey weight composition +# =========================================================================== + +cat("\n--- S2: StackedDiD (stacking + Q-weight + composed weights) ---\n") + +kappa_pre <- 1L +kappa_post <- 1L +cohorts <- c(3, 4) # adoption events + +# Step 1: Build sub-experiments +sub_experiments <- list() +for (a in cohorts) { + # Treated units: first_treat == a + treated_units <- unique(staggered_data$unit[staggered_data$first_treat == a]) + # Control units: never-treated + control_units <- unique(staggered_data$unit[is.infinite(staggered_data$first_treat)]) + + if (length(treated_units) == 0 || length(control_units) == 0) next + + # Time window: [a - kappa_pre, a + kappa_post] + t_start <- a - kappa_pre + t_end <- a + kappa_post + + all_units <- c(treated_units, control_units) + sub_df <- staggered_data[staggered_data$unit %in% all_units & + staggered_data$period >= t_start & + staggered_data$period <= t_end, ] + + if (nrow(sub_df) == 0) next + + sub_df$sub_exp <- a + sub_df$event_time <- sub_df$period - a + sub_df$D_sa <- as.integer(sub_df$unit %in% treated_units) + + sub_experiments[[as.character(a)]] <- sub_df +} + +stacked <- do.call(rbind, sub_experiments) +rownames(stacked) <- NULL +n_stacked <- nrow(stacked) +cat(sprintf(" Stacked dataset: %d observations\n", n_stacked)) + +# Step 2: Compute Q-weights (sample_share) +# Count distinct units per sub-experiment +N_D <- sapply(sub_experiments, function(df) length(unique(df$unit[df$D_sa == 1]))) +N_C <- sapply(sub_experiments, function(df) length(unique(df$unit[df$D_sa == 0]))) +N_Omega_D <- sum(N_D) +N_Omega_C <- sum(N_C) +N_grand <- N_Omega_D + N_Omega_C + +cat(sprintf(" N_D per sub-exp: %s\n", paste(N_D, collapse = ", "))) +cat(sprintf(" N_C per sub-exp: %s\n", paste(N_C, collapse = ", "))) +cat(sprintf(" N_Omega_D = %d, N_Omega_C = %d, N_grand = %d\n", + N_Omega_D, N_Omega_C, N_grand)) + +q_control <- setNames(numeric(length(cohorts)), as.character(cohorts)) +for (a_str in names(N_D)) { + n_c <- N_C[a_str] + n_d <- N_D[a_str] + if (n_c == 0 || N_Omega_C == 0) { + q_control[a_str] <- 1.0 + } else { + control_share <- n_c / N_Omega_C + sample_share <- (n_d + n_c) / N_grand + q_control[a_str] <- sample_share / control_share + } +} +cat(sprintf(" Q-weights (control): %s\n", + paste(sprintf("%s=%.6f", names(q_control), q_control), collapse = ", "))) + +# Assign Q-weights: treated=1, control=q_control[sub_exp] +stacked$Q_weight <- 1.0 +for (a_str in names(q_control)) { + mask <- stacked$sub_exp == as.numeric(a_str) & stacked$D_sa == 0 + stacked$Q_weight[mask] <- q_control[a_str] +} + +# Step 3: Compose Q x survey weight, normalize +stacked$composed_weight <- stacked$Q_weight * stacked$weight +w_sum <- sum(stacked$composed_weight) +if (w_sum > 0) { + stacked$composed_weight <- stacked$composed_weight * (n_stacked / w_sum) +} + +# Step 4: Fit svyglm +# Use strata/PSU structure from original data (propagated into stacked data) +# to match Python's TSL variance on the re-resolved stacked survey design. +# Omit FPC — the stacked data inflates the sample size per stratum due to +# unit duplication across sub-experiments, which can make n_h > N_h. +stacked$event_time_fac <- relevel(factor(stacked$event_time), ref = "-1") +design_s2 <- svydesign(ids = ~psu, strata = ~stratum, + weights = ~composed_weight, data = stacked) + +fit_s2 <- svyglm(outcome ~ D_sa + event_time_fac + D_sa:event_time_fac, + design = design_s2) + +coefs_s2 <- coef(fit_s2) +ses_s2 <- SE(fit_s2) +V_s2 <- vcov(fit_s2) +df_s2 <- degf(design_s2) + +# Extract event-study interaction coefficients (D_sa:event_time_fac*) +interaction_names <- grep("^D_sa:event_time_fac", names(coefs_s2), value = TRUE) +event_study_s2 <- list() +post_coef_indices <- integer(0) + +for (nm in interaction_names) { + # Parse event time from name like "D_sa:event_time_fac0" or "D_sa:event_time_fac1" + e <- as.integer(sub("^D_sa:event_time_fac", "", nm)) + event_study_s2[[as.character(e)]] <- list( + coef = as.numeric(coefs_s2[nm]), + se = as.numeric(ses_s2[nm]) + ) + if (e >= 0) { + post_coef_indices <- c(post_coef_indices, which(names(coefs_s2) == nm)) + } +} + +# ATT = mean of post-treatment interaction coefficients +post_coefs <- coefs_s2[post_coef_indices] +n_post <- length(post_coefs) +att_s2 <- mean(post_coefs) + +# SE via delta method: ATT = (1/n_post) * sum(post_coefs) +# var(ATT) = (1/n_post)^2 * ones' V_sub ones +V_sub_s2 <- V_s2[post_coef_indices, post_coef_indices] +ones <- rep(1, n_post) +se_s2 <- sqrt(as.numeric(t(ones) %*% V_sub_s2 %*% ones)) / n_post + +results$s2_stacked_did <- list( + scenario = "S2", + description = "StackedDiD: stacked WLS with Q-weight x survey weight composition", + estimator = "survey::svyglm", + att = att_s2, + se = se_s2, + n_stacked = n_stacked, + df = df_s2, + event_study = event_study_s2 +) + +cat(sprintf(" ATT = %.6f (SE = %.6f)\n", att_s2, se_s2)) +cat(sprintf(" Event-study effects: %s\n", + paste(sapply(names(event_study_s2), function(e) + sprintf("e=%s: %.4f", e, event_study_s2[[e]]$coef)), collapse = ", "))) + +# =========================================================================== +# S3: SunAbraham — interaction-weighted ATT +# =========================================================================== + +cat("\n--- S3: SunAbraham (IW-aggregated ATT) ---\n") + +# Create cohort x relative-time interaction dummies +# Cohort 3: rel_times = {-2, -1, 0, 1, 2} -> exclude -1 -> {-2, 0, 1, 2} +# Cohort 4: rel_times = {-3, -2, -1, 0, 1} -> exclude -1 -> {-3, -2, 0, 1} +staggered_data$rel_time_3 <- ifelse(staggered_data$first_treat == 3, + staggered_data$period - 3, NA) +staggered_data$rel_time_4 <- ifelse(staggered_data$first_treat == 4, + staggered_data$period - 4, NA) + +# Build interaction dummies +interaction_info <- list( + list(g = 3, e = -2, name = "D_3_m2"), + list(g = 3, e = 0, name = "D_3_0"), + list(g = 3, e = 1, name = "D_3_1"), + list(g = 3, e = 2, name = "D_3_2"), + list(g = 4, e = -3, name = "D_4_m3"), + list(g = 4, e = -2, name = "D_4_m2"), + list(g = 4, e = 0, name = "D_4_0"), + list(g = 4, e = 1, name = "D_4_1") +) + +for (info in interaction_info) { + staggered_data[[info$name]] <- as.integer( + staggered_data$first_treat == info$g & + (staggered_data$period - info$g) == info$e + ) +} + +# Verify non-zero support for each interaction +for (info in interaction_info) { + n_ones <- sum(staggered_data[[info$name]]) + cat(sprintf(" %s: %d observations\n", info$name, n_ones)) +} + +design_s3 <- svydesign(ids = ~psu, strata = ~stratum, fpc = ~fpc, + weights = ~weight, data = staggered_data) + +# Formula: outcome ~ factor(unit) + factor(period) + all interaction dummies +interaction_names_s3 <- sapply(interaction_info, function(x) x$name) +formula_rhs <- paste(c("factor(unit)", "factor(period)", interaction_names_s3), + collapse = " + ") +formula_s3 <- as.formula(paste("outcome ~", formula_rhs)) + +fit_s3 <- svyglm(formula_s3, design = design_s3) + +coefs_s3 <- coef(fit_s3) +ses_s3 <- SE(fit_s3) +V_s3 <- vcov(fit_s3) +df_s3 <- degf(design_s3) + +# Extract cohort effects +cohort_effects_s3 <- list() +for (info in interaction_info) { + nm <- info$name + cohort_effects_s3[[nm]] <- list( + coef = as.numeric(coefs_s3[nm]), + se = as.numeric(ses_s3[nm]) + ) + cat(sprintf(" delta(%d, %d) = %.6f (SE = %.6f)\n", + info$g, info$e, coefs_s3[nm], ses_s3[nm])) +} + +# --- IW Aggregation --- +# Post-treatment interactions (e >= 0) +post_info <- Filter(function(x) x$e >= 0, interaction_info) + +# Compute survey-weighted mass n_{g,e} for each post-treatment (g, e) pair +post_masses <- list() +for (info in post_info) { + mask <- staggered_data$first_treat == info$g & + (staggered_data$period - info$g) == info$e + n_ge <- sum(staggered_data$weight[mask]) + post_masses[[info$name]] <- list(g = info$g, e = info$e, mass = n_ge, + coef_name = info$name) +} + +# Get unique post-treatment event times +post_e_vals <- sort(unique(sapply(post_info, function(x) x$e))) + +# Per-period IW aggregation: beta_e = sum_g w_{g,e} * delta_{g,e} +period_effects <- list() +for (e in post_e_vals) { + pairs_at_e <- Filter(function(x) x$e == e, post_masses) + masses_at_e <- sapply(pairs_at_e, function(x) x$mass) + total_mass_at_e <- sum(masses_at_e) + + if (total_mass_at_e == 0) next + + beta_e <- 0.0 + for (p in pairs_at_e) { + w_ge <- p$mass / total_mass_at_e + delta_ge <- as.numeric(coefs_s3[p$coef_name]) + beta_e <- beta_e + w_ge * delta_ge + } + period_effects[[as.character(e)]] <- list(effect = beta_e, mass = total_mass_at_e) +} + +# Overall ATT = sum_e w_e * beta_e, where w_e proportional to total mass at e +total_post_mass <- sum(sapply(period_effects, function(x) x$mass)) +att_s3 <- 0.0 +for (e_str in names(period_effects)) { + w_e <- period_effects[[e_str]]$mass / total_post_mass + att_s3 <- att_s3 + w_e * period_effects[[e_str]]$effect +} + +# SE via delta method: ATT = w' * delta +# Build the full weight vector for each (g, e) in post-treatment +overall_weights <- numeric(length(post_info)) +overall_coef_indices <- integer(length(post_info)) + +for (i in seq_along(post_info)) { + info <- post_info[[i]] + nm <- info$name + e <- info$e + mass_ge <- post_masses[[nm]]$mass + + # Per-period mass + e_str <- as.character(e) + total_at_e <- period_effects[[e_str]]$mass + w_ge_within_e <- mass_ge / total_at_e # cohort weight within period + + # Period weight + w_e <- period_effects[[e_str]]$mass / total_post_mass + + overall_weights[i] <- w_e * w_ge_within_e + overall_coef_indices[i] <- which(names(coefs_s3) == nm) +} + +V_sub_s3 <- V_s3[overall_coef_indices, overall_coef_indices] +var_att_s3 <- as.numeric(t(overall_weights) %*% V_sub_s3 %*% overall_weights) +se_s3 <- sqrt(max(var_att_s3, 0)) + +cat(sprintf(" IW-aggregated ATT = %.6f (SE = %.6f)\n", att_s3, se_s3)) + +# Output with g_e keys (negative as m) +cohort_effects_output <- list() +for (info in post_info) { + e_label <- if (info$e < 0) paste0("m", abs(info$e)) else as.character(info$e) + key <- paste0("g", info$g, "_e", e_label) + cohort_effects_output[[key]] <- cohort_effects_s3[[info$name]] +} + +results$s3_sun_abraham <- list( + scenario = "S3", + description = "SunAbraham: IW-aggregated ATT with survey-weighted cohort masses", + estimator = "survey::svyglm", + att = att_s3, + se = se_s3, + df = df_s3, + cohort_effects = cohort_effects_output +) + +# =========================================================================== +# S4: TripleDifference — three-way interaction +# =========================================================================== + +cat("\n--- S4: TripleDifference (three-way interaction DDD) ---\n") + +design_s4 <- svydesign(ids = ~psu, strata = ~stratum, fpc = ~fpc, + weights = ~weight, data = ddd_data) + +fit_s4 <- svyglm(outcome ~ group * partition * time, design = design_s4) + +s4_res <- extract_svyglm_results(fit_s4, design_s4, "group:partition:time") + +results$s4_triple_diff <- list( + scenario = "S4", + description = "TripleDifference: three-way interaction DDD coefficient", + estimator = "survey::svyglm", + att = s4_res$att, + se = s4_res$se, + t_stat = s4_res$t_stat, + df = s4_res$df, + ci_lower = s4_res$ci_lower, + ci_upper = s4_res$ci_upper +) + +cat(sprintf(" DDD = %.6f (SE = %.6f)\n", s4_res$att, s4_res$se)) +cat(sprintf(" df = %d, CI = [%.4f, %.4f]\n", s4_res$df, s4_res$ci_lower, s4_res$ci_upper)) + +# =========================================================================== +# Embed data and metadata +# =========================================================================== + +# Staggered data: export all columns as lists +staggered_export <- as.list(staggered_data[, c("unit", "period", "stratum", "psu", "fpc", + "weight", "first_treat", "x1", "x2", + "outcome")]) + +# DDD data: export all columns as lists +ddd_export <- as.list(ddd_data[, c("obs_id", "group", "partition", "time", + "outcome", "stratum", "psu", "fpc", "weight")]) + +results[["_staggered_data"]] <- staggered_export +results[["_ddd_data"]] <- ddd_export +results[["_metadata"]] <- list( + description = "Golden values for survey-enabled estimator validation", + r_version = paste0(R.version$major, ".", R.version$minor), + survey_version = as.character(packageVersion("survey")), + seed = 42, + n_units = length(unique(staggered_data$unit)), + n_periods = length(unique(staggered_data$period)), + n_obs_staggered = nrow(staggered_data), + n_obs_ddd = nrow(ddd_data) +) + +# =========================================================================== +# Save +# =========================================================================== + +json_path <- file.path(output_dir, "survey_estimator_validation_golden.json") +writeLines(toJSON(results, auto_unbox = TRUE, pretty = TRUE, digits = 10), json_path) +cat(sprintf("\nResults saved to: %s\n", json_path)) + +cat("\n=== Summary ===\n") +cat(sprintf("S1 (ImputationDiD): coef_x1 = %.6f, SE = %.6f\n", + results$s1_imputation_did$coef_x1, results$s1_imputation_did$se_x1)) +cat(sprintf("S2 (StackedDiD): ATT = %.6f, SE = %.6f\n", + results$s2_stacked_did$att, results$s2_stacked_did$se)) +cat(sprintf("S3 (SunAbraham): ATT = %.6f, SE = %.6f\n", + results$s3_sun_abraham$att, results$s3_sun_abraham$se)) +cat(sprintf("S4 (TripleDifference): DDD = %.6f, SE = %.6f\n", + results$s4_triple_diff$att, results$s4_triple_diff$se)) +cat("Done.\n") diff --git a/benchmarks/data/synthetic/survey_estimator_validation_golden.json b/benchmarks/data/synthetic/survey_estimator_validation_golden.json new file mode 100644 index 00000000..bf8a6d71 --- /dev/null +++ b/benchmarks/data/synthetic/survey_estimator_validation_golden.json @@ -0,0 +1,106 @@ +{ + "s1_imputation_did": { + "scenario": "S1", + "description": "ImputationDiD: control-only WLS regression on Omega_0", + "estimator": "survey::svyglm", + "coef_x1": 2.030610978, + "se_x1": 0.037894492356, + "coef_x2": 1.397412961, + "se_x2": 0.13649612979, + "df": 6, + "n_obs_omega0": 558 + }, + "s2_stacked_did": { + "scenario": "S2", + "description": "StackedDiD: stacked WLS with Q-weight x survey weight composition", + "estimator": "survey::svyglm", + "att": 5.039606586, + "se": 0.43170156311, + "n_stacked": 678, + "df": 6, + "event_study": { + "0": { + "coef": 5.0048694773, + "se": 0.4629603277 + }, + "1": { + "coef": 5.0743436947, + "se": 0.56222065855 + } + } + }, + "s3_sun_abraham": { + "scenario": "S3", + "description": "SunAbraham: IW-aggregated ATT with survey-weighted cohort masses", + "estimator": "survey::svyglm", + "att": 5.2043029801, + "se": 0.35214822631, + "df": 6, + "cohort_effects": { + "g3_e0": { + "coef": 4.8522956607, + "se": 0.50822504026 + }, + "g3_e1": { + "coef": 4.450329809, + "se": 0.80872622454 + }, + "g3_e2": { + "coef": 5.8074158445, + "se": 0.34262341419 + }, + "g4_e0": { + "coef": 5.2165186904, + "se": 0.51731434905 + }, + "g4_e1": { + "coef": 5.93122286, + "se": 0.55519324739 + } + } + }, + "s4_triple_diff": { + "scenario": "S4", + "description": "TripleDifference: three-way interaction DDD coefficient", + "estimator": "survey::svyglm", + "att": 2.9254016361, + "se": 0.61785760375, + "t_stat": 4.7347505611, + "df": 6, + "ci_lower": 1.4135585431, + "ci_upper": 4.437244729 + }, + "_staggered_data": { + "unit": [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150], + "period": [1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5], + "stratum": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], + "psu": [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9], + "fpc": [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250], + "weight": [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 4.4444444444, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5], + "first_treat": [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf", "Inf"], + "x1": [0.75930152485, -0.45154482235, 2.5671695492, 1.4263542986, 1.5321546753, -0.43825232799, 1.0922379662, 0.59971782205, 1.16650767, 3.2876752456, 0.62447412976, 0.028217615849, 1.5632068825, 1.3372615844, 2.590013548, 2.2324219277, 0.10730748889, 1.8546502976, 1.7848980587, 2.8955655823, 0.57022186105, 1.9366190685, -0.34505586134, 3.3894589213, 2.8647779785, 0.64922401111, -0.34900713014, 2.0430087, 2.1831938638, 1.9936443736, 0.95145589286, 0.51589102965, 1.7688067326, 1.9946543397, 1.7207406267, -0.53623665489, 1.8007488184, 1.9541966223, 0.8636934072, 0.40541183799, 1.0049585806, 0.75491202203, 1.6526117034, 0.40599753452, 1.0408295556, 1.8857748537, 1.5037749047, 1.9864875367, 3.5152284462, 2.1288214286, -1.2009292377, 1.4337771974, 2.5713251274, 3.7595392423, 0.62313840176, -0.35085556563, 0.39417860524, 0.34594421792, 1.0542562769, 1.8146220323, -0.90122205074, 2.636972167, 1.0077747449, 1.1158918995, 1.9956196416, 0.33741518612, 0.46791196304, 2.3767874236, 0.98296978991, 0.21639779591, 0.68566789044, 0.24848712647, 0.37820390664, 0.13186879931, 1.9283659033, 0.12598176557, 1.1156677286, 0.66563472269, 0.54149657418, 2.7502366041, 0.028236284889, 1.5479519959, -0.30158243011, 0.73388390362, 1.2306486048, -0.090965408131, 1.948707012, 0.8772352987, 1.4442258511, 0.55762829214, -0.42921727651, 1.5980689086, 2.1584816646, 2.4488636888, 0.11936295047, 2.3499606936, 1.6168728298, 0.87328253586, 1.9036077788, 0.52861477085, -0.79615624156, 0.64905045094, -0.29849585656, 1.3900189986, 2.7977058996, -0.73387372296, -0.13844075421, 0.94656393945, 0.18240388017, 1.1167160401, 1.1727554117, 1.569545014, 1.283846665, -0.65155566307, 1.4460032632, 1.3647732143, 1.4131950374, 0.7091835259, -1.4999298086, 1.5609666388, 0.87375169748, 0.64580357973, 1.0574125402, 1.6315653729, 1.1034502639, 2.1099782258, 1.5703933999, 0.15732972943, 3.0815754564, 3.204458937, 1.6240739637, -0.56262940219, 0.83069365639, 2.3355138173, 2.0437220076, 1.1480123037, 3.5595935489, 0.5816196756, -0.41320011493, 2.2736952724, 0.11240315876, 1.546041053, 0.58761527621, 3.9120554803, 1.8762940284, 0.32266449397, 0.93373850851, 2.2625633836, 1.7973404852, 0.37438326079, 0.79537923217, 1.8602421677, 1.4389909129, 2.4350721416, 1.853527373, 0.74211266462, 1.5823694661, 2.3929436368, 0.45360450197, 1.9665124752, 0.72903781875, 0.34107934627, 0.36564063906, 1.069268046, 2.58680772, 0.38367734382, 0.31511219049, 1.5634163195, 0.46328576489, 3.0458737762, 0.31540458378, 1.2891288117, 1.4510063316, 1.6997593311, 3.8093820424, -0.02532795711, 2.2454704522, 1.4315731876, 0.86479419469, 1.931236351, 1.546771687, 0.67448126561, 0.62791776481, 1.8527641067, 2.9885968452, 0.72654166533, -0.28702655356, 0.093324095625, 0.93160467492, 1.4728918746, 0.77857293496, 1.770498071, 0.96538296143, 0.58612021669, 2.6071059949, 1.0754569687, 2.2573470694, -0.28248085948, 1.7873190889, 3.3533618939, 1.0241738007, -0.23255282576, 1.6325284868, 0.32807313, 1.0466024888, 1.4875342786, 0.30985468821, 1.728546145, 0.90877229122, -0.076362404748, -0.54881569695, -0.48851986196, 0.41570942888, 0.86368879099, 1.3466421093, 0.056752771425, 2.4922020416, -0.48599833729, 0.78517569929, 1.8490815281, 1.9284422658, 0.68852189571, 2.1391507084, -0.44455553575, 2.9463565254, -0.39056017167, 0.32356891455, -0.2094187599, 1.3338693164, 3.2853390517, 2.7221633553, -0.47682890212, 1.3859411104, 2.5885217387, 1.3043431827, 0.081825202294, 0.29522204545, 1.4978327241, 2.5974294108, 2.1876197612, 0.6201880143, 0.29813007422, 1.3983486856, 0.65046308208, 1.2207434964, 1.3965134431, 1.0420130882, 1.141016294, 0.94439234464, 2.4310329015, 1.6349125854, -0.26927176388, 0.95548695467, 1.2490669132, 0.92164427164, -0.698738656, 0.59756721999, 1.5555118828, 2.676842279, -0.40915278828, -0.40243947331, 0.28856978199, -0.76315703067, 0.44946655821, 0.72264824082, -0.42256969986, -1.5888345991, 1.1134185503, 0.56807706375, 3.0204911938, 1.0959559494, -0.85352956513, 0.99839542101, 0.60622901577, 2.3882811686, 0.8530704155, 0.54297637416, 1.8383970362, 1.8526081586, 1.8353824183, 2.8198906215, 0.57061411368, 0.92921302702, 0.15406307598, 1.9594732767, 1.6903563055, -0.97138785093, 1.1499348799, 0.51834957335, 3.441937265, 2.1578955389, 1.4345028473, 2.8293380804, 0.83268214912, 2.950651725, 0.21498849137, 1.4209575225, 1.1006039833, 1.4214569167, 2.5461151583, -0.50382115364, 0.84908553496, 1.5710073737, 1.2965325213, 2.1046594408, 0.48111921327, 2.7183439359, 2.1141886015, 4.6658653698, 1.2049223954, 1.6143659147, 3.1980308101, 1.7009800552, 0.61692485816, 0.9936774982, 0.76458543508, 2.4091243609, 2.1504004866, -0.43836832761, 1.2996458911, 0.79094352471, -0.35813348724, 2.0945296465, -0.76133547519, 2.143289764, 0.40877788087, 0.60883591439, 1.1163525475, 2.0147948047, 2.3963265779, 0.57439628896, -0.82495043034, -0.039229298322, 0.23034222603, 2.7618634469, 0.55638501813, 1.3696766074, -0.15892752096, 1.1644119933, 2.5624519732, 0.62167387658, 0.9848640138, 1.3279385279, 2.9109098069, 1.3851030992, 1.4761264583, 1.9985996055, 0.21068209648, 1.8212588498, 1.9887783139, 1.8291407193, 2.0147748677, 1.3975437333, 1.8360095524, 1.2798464551, 0.60187566995, 0.070791193952, 0.43304410415, 0.47918691108, 2.8362077039, 1.9149717985, 0.68705514707, 0.27102260199, 1.6120686682, 4.2419036933, 2.8413131708, -0.61980033672, 1.0430933439, 3.2332804922, 1.9617559034, 2.3974132423, 0.76641460669, 2.0049859239, 1.9242412393, 5.2290694955, 1.7204525668, -0.10653906922, 0.79611883084, 2.0702352846, 0.09899935332, -1.004409945, -0.021822606355, 1.0520703934, 2.9389021486, 1.7258024351, 0.46238071556, 0.53539309025, 0.19476290297, 2.5621972688, 0.40348681192, 0.071566481453, 0.25217191939, 1.4321286195, 2.8072345623, 2.0138145261, 1.6823731606, 1.3630970485, 0.27541545537, 1.2815438004, 2.8760791114, -1.2612657862, 0.92477066059, 0.74320968305, 2.0777838607, 2.2501665752, 0.60105340384, 2.0928111172, -0.62549380017, 2.1107175585, -0.079532027858, 0.88771625712, 0.689642296, 1.8672616569, 1.2788123529, -0.068699836486, -1.707823184, 1.1409735242, 0.82662340189, 0.62898160644, 1.1889315401, -0.3714150671, 0.44265954792, -0.031305071, -0.78300947947, 1.2804001958, 1.3452758695, 2.4773295663, 0.90260619596, 1.1193300653, 2.4629822812, 0.35357101655, 0.16510159109, -0.83729725534, -0.063695599782, 1.9063085121, -1.1596539678, 1.6484573702, -0.44643053945, 1.0064294419, 1.4976640429, 0.28717026613, 0.75194732198, 1.4985111314, 1.0737875627, 1.2514644342, 0.4603273949, 0.1663580584, 2.4374124191, -0.97024657662, 2.5270046193, 0.05151706681, 1.0223203858, 1.8876532936, 2.035568172, 0.83947814095, 1.8640694933, -1.0229759353, 1.7638963735, 0.68839722656, -0.41736502545, -1.5658138486, 0.84517906344, 3.1235343019, 1.4733760641, 2.6307847945, 1.1386693741, -0.054615277878, 1.8539609785, 1.5529514655, 1.7065992926, 1.3011132343, 1.3474519916, 0.27342516328, 1.5952236859, 0.60783204595, 0.93081841009, 0.16729481672, 1.3521386442, 1.5679990452, 1.2296125103, 0.76551261771, 1.1743562256, 0.66930000984, 2.3722425083, 2.8927027067, -0.36199125571, -0.17736920673, 1.4135385789, 0.28668776196, 1.0505761947, 1.1029326991, 0.026523149019, -1.0281251678, 1.8643908338, -0.10254022379, -1.0546002575, -2.4179326794, 1.7312378222, 1.4510970888, 1.9622934659, 1.6447922232, 1.0580284757, 0.2944240943, 2.2637756523, 3.3033647558, -0.70022093779, 0.49301076488, 1.1077549343, 0.41031665668, 2.6941058488, 0.2008180971, 2.356906644, 1.4535080136, 2.4280925177, 3.5610980682, 1.0656247536, 2.1767262592, 1.6106979317, 0.18832647418, 2.0224022599, 1.518136206, 1.5894570178, 1.2261116541, 0.48230051133, 2.6463983687, -0.11645602835, -0.15182344243, 1.9949277184, 0.46718933481, 2.2443644149, 0.80277218534, -0.22820968603, 2.5796964118, 1.1071950009, 3.1999783163, 0.32496632039, 0.52494290385, 1.3687739754, 1.341943004, 1.6433989224, -0.077664383434, -0.11289712518, 2.0132866187, 0.89379665854, 0.62354307009, 0.29215210084, 0.29906451326, -0.79278568567, 1.4090628507, 2.1671741212, 1.0946923565, 1.4927412654, 0.399156287, 1.3742728802, 0.99165119457, 0.16456851797, -0.10984068799, 0.28353619935, 2.3298811626, 1.7274784298, 0.54115883162, 2.8295581799, 1.3416078346, 1.1629362154, 2.7472866962, 0.31274216461, 2.4729077815, 1.0223276394, 1.0838473859, 0.83187494893, 1.1286403586, 2.5665106274, 1.0439904549, 1.9614676419, 2.3333288552, 2.2221932403, 1.7638766007, 0.3263448443, 1.0030982253, 1.2538695431, 0.94157287213, 1.0960527793, 1.7679375011, 1.0426570853, 1.6236533298, 1.5413600604, 1.0003932436, 0.74571012139, 2.6711643741, 2.0134963116, -0.11653466702, 2.8096885613, 0.23189901005, -0.081036260153, -0.25313244857, 1.4511259695, 0.5671673744, 1.124440643, 1.9896269726, 1.5335159382, -0.80806017494, -0.84978379604, 1.0593000818, 1.8747255273, -0.27777756078, 1.0905244809, 1.5223063575, 2.6947373632, 1.8647140186, 2.2049538408, 1.4046036267, 1.0808592308, 1.3202856526, 1.8159844644, 2.744173364, 0.36870939647, 0.60047104809, 0.53483792324, 0.74224271673, 2.3267999713, 2.3473722257, 0.13114043114, 1.2115595737, 0.66999883854, 2.9080864419, -0.017381864708, 0.77900396654, 1.2480281877, 0.14572097873, 1.3952557119, 0.071657382905, 1.2753555624, -0.33324469351, 3.8455577221e-05, 2.2658665766, -0.28809786128, -0.060295829976, 1.0130135222, 0.87960120031, 3.3663814969, 0.55953145577, 0.76156007092, 1.8310749093, 1.1400532507, 1.5487397001, -0.77287540073, -1.6929714305, -0.30720685006, 1.3141094298, 0.46670292089, 0.98880775234, 1.3250830177, 1.1173802109, 0.99834326838, 0.13431013941, -0.0089376086772, 0.14709711346, 1.5632291445, 2.5086295367, 2.0010403099, -0.8282885352, 2.2709973054, 1.9103530325, 1.4235212159, -0.70648459321, -0.65458561855, 0.53142693279, 1.6613059016, 0.020095800042, 4.7111989568, -2.253824851, 0.36406605067, 0.64043735674, 0.53663308338, 1.1810092892, 1.0423951576, -0.27429267579, -1.1828137663, 1.2937679028, 1.498180188, 0.28689934763, 1.2679680746, 0.88683229705, 1.9760470305, -0.51073533264, -0.82818049497, 0.94879957159, 0.5471015179, 2.1447745979, 0.49528014747, 1.0239027325, -0.06883274099, -0.21304019372, 0.85719481063, 1.5497793232, -0.92768196153, -0.16400591083, -0.34618206162, 2.2167739778, 2.2236126241, -0.73252713711, 1.1573460041, 0.6444185769, 0.086608242853, 1.9211974396, -0.022950311846, 1.4801799012, 0.70513639322, 2.3880549821, 0.99023548823, 0.51954801087, 0.97029442136, 1.1791083937, 1.4019424373, 1.4870034576, 0.20913499392, 1.9650313653, 1.8081305689, 0.5901118286, 2.8964372212, 0.44479837331, -0.040606771595, 1.0698022596, 1.0428125126, 1.6009395945, -0.67339415415, -0.21982527525, 2.262918969, 2.1613709305, 0.6162755256, -0.60009211544, 2.3233331328, 2.8090421621, 0.42285880414, 0.19769473739, 2.9234951668, 0.82962891435, 1.0867498921, 1.276136635, 2.9048601863, 0.10827776466, 2.0593921642, 0.67934454887, 1.7051458863, 0.46650299393, 0.47047347146, 1.8006681609, 0.73659408523, 2.4824758631, 4.2271963879], + "x2": [0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0], + "outcome": [19.4739316745, 23.3641631197, 40.7412791806, 40.8857192998, 48.1034439461, 13.0350853756, 21.73675554, 31.8596799445, 36.5434061086, 44.6285614176, 18.446455767, 24.2834737055, 32.7234045243, 39.8729585835, 47.1792150022, 22.2128796616, 20.6372870132, 36.4088304776, 41.6467548256, 46.9442661461, 17.6077192073, 24.5149132831, 31.3272983888, 44.2578739228, 47.2305122396, 16.4006520535, 18.4313474168, 36.4632763167, 39.6629246366, 43.6391330894, 21.3171730436, 24.8097262146, 37.8368461385, 43.5873078453, 46.318570446, 15.5122179785, 24.178147476, 29.5513630107, 36.1326112151, 40.600668923, 20.6378306765, 24.1702851652, 34.5695012089, 41.082973278, 47.0652821331, 18.4169313313, 24.5935836002, 27.197738972, 42.6492494796, 43.8181933351, 15.0925699916, 24.8663193605, 35.3490849422, 44.001650182, 46.3186256044, 20.2780436704, 25.219998487, 29.1382750287, 42.2733488196, 47.2462257305, 9.7283046099, 22.0508701644, 23.8933438592, 31.657760582, 35.6105121162, 15.206603251, 22.6157410297, 27.6556178689, 30.7187297535, 35.7350115184, 17.2281733674, 17.6952481098, 24.6319559119, 30.8131213328, 38.1635002245, 17.6192274872, 24.8574670292, 29.8226366096, 33.3638031923, 43.7483324427, 16.561120395, 21.5085385712, 25.2664570625, 31.8609568628, 37.309746093, 7.9968816858, 20.1236731099, 24.2875778442, 28.2326850829, 30.1324715147, 8.3521441417, 18.1786633182, 25.4349915087, 31.6294319673, 29.9825744807, 23.9986947932, 25.8895350223, 29.2861883019, 38.2234153382, 40.252408873, 12.563836974, 20.6716074231, 24.1746802947, 32.513484526, 41.0505218505, 11.1198137254, 17.6239694307, 23.1947095892, 29.246808938, 34.6628565739, 16.3010063569, 23.3041966446, 26.5182405659, 28.677278069, 35.7669748631, 21.9695033489, 25.1106291122, 29.4902662038, 29.4278304189, 42.460445638, 20.0418321702, 24.1218861231, 31.6833442275, 39.3146382846, 40.7635752665, 19.2363128401, 22.0064770976, 29.4080338713, 39.9073403956, 46.9199233444, 18.8298468937, 18.6049845695, 30.8813648526, 41.7097913591, 43.5945989437, 14.6431108687, 25.051808905, 27.8368664006, 30.1139679992, 42.2550320568, 18.7291488413, 25.331531336, 33.8469805815, 45.6138479055, 45.0798170259, 16.3578442553, 20.6421198442, 34.0364629392, 37.256420466, 40.0053549474, 17.5950520074, 24.2834834483, 33.6705476341, 39.6819535725, 44.3348785801, 18.3041653308, 25.6687480626, 36.7224692076, 39.0573036602, 47.741062082, 19.0162508084, 24.3643571017, 35.2998016692, 40.4004599448, 45.8444349168, 16.9416496735, 21.6903306005, 32.8508944045, 35.2999160667, 48.7828256461, 14.0464692399, 25.0950309075, 29.0300295787, 40.7597315824, 49.3765476193, 11.6393615345, 20.2788687991, 25.8391623225, 33.001659223, 41.0391339888, 15.8581152774, 20.2201616362, 25.5810788137, 37.5797687759, 46.4413821463, 14.6991755277, 16.7456422392, 23.5120513581, 33.4458218445, 40.7390805344, 11.4238078891, 20.0586590775, 23.3448139322, 30.8205378328, 39.7225788909, 19.170380191, 25.7674883707, 25.4635046449, 39.8672099491, 47.844201602, 19.0281735509, 20.8037794699, 28.6975759224, 30.743959002, 39.3625982001, 17.8437179552, 18.1169680619, 27.3008689933, 30.3026802736, 35.532195355, 17.6226595308, 21.1447971262, 26.0652385371, 34.4994882275, 39.7449807105, 14.5678797698, 25.2450680427, 23.455616794, 31.1194867427, 38.3846639066, 18.8563582951, 18.9358832127, 28.2135003516, 26.918937557, 36.9457418409, 14.4831327106, 22.8769577443, 25.9595257471, 34.3652996993, 41.6412546388, 20.6620049701, 19.7939856637, 26.4669348632, 36.0771511009, 35.3774654368, 19.7816155657, 24.9903574059, 32.9393059608, 40.4979401463, 43.2100454216, 15.6758597129, 22.3397472829, 27.1046553811, 31.4397278878, 35.1692907137, 18.6055416283, 26.0180787957, 30.6039939209, 34.6092372619, 43.652127249, 18.6698670443, 20.7271159619, 28.7102292199, 30.786759811, 37.959526707, 11.5233646115, 20.5538543585, 26.5471430876, 36.6285576057, 34.4094780532, 18.1883630366, 25.8576589018, 28.2267393777, 36.2215090307, 40.9584351467, 14.1913469586, 18.4557646574, 27.7348220528, 32.9357435577, 42.2672396555, 18.7106737016, 18.0521006426, 28.1649031515, 31.8223655336, 38.8424401635, 19.5396942232, 23.1841593164, 34.121406408, 40.3670935018, 44.0368668268, 22.5456544199, 24.5072199716, 34.7358197851, 35.2477934105, 47.3756961169, 20.8467580058, 19.6022570254, 33.9631490154, 37.6255053281, 48.958379002, 14.5494860395, 16.8377483155, 32.6519838393, 32.1035003452, 42.3328717408, 15.6549641761, 23.8346782961, 35.2725786011, 40.7086837812, 46.8361181083, 14.7267445358, 22.2267928513, 33.3310385798, 37.3694039503, 45.2504958745, 17.0690384627, 24.4550370173, 35.2538627429, 46.3118386298, 43.4917989547, 18.5085031716, 26.4871944632, 34.1074665836, 38.0631441866, 43.8632306635, 19.7715987522, 29.8773096884, 38.1313969619, 36.7500270215, 47.5808153362, 16.3457256022, 18.3675263326, 32.8568250293, 34.7203315673, 42.6322803536, 18.4519576786, 27.257071429, 37.1753012344, 42.9493285424, 49.4310675977, 20.1054450594, 20.7114231987, 31.5241501307, 37.5842736361, 49.5644165859, 17.2052273399, 25.5351221761, 31.8885044524, 41.6074774929, 49.2286139756, 18.0642917908, 24.2494362115, 28.2868056267, 42.7588568066, 46.8324124922, 19.4800447456, 25.3073384954, 25.1870685398, 41.0147760346, 47.1115181044, 16.0857206407, 22.0245958713, 27.5937468363, 37.9482792567, 42.4010481422, 15.7803141181, 21.1428567868, 26.9316090416, 37.9168884791, 46.3131080178, 22.4124590076, 21.6378390605, 27.581624088, 40.1153227618, 52.8685090406, 21.5022278006, 17.8856469927, 24.5009089085, 42.7559165401, 45.0593119413, 20.4236964476, 23.64785988, 29.416327389, 38.4449244846, 50.8218711915, 21.1828192155, 20.889343759, 29.0827327003, 39.8838564718, 43.35127852, 15.4963091545, 21.9749323807, 30.6767893026, 42.1348375651, 46.4444550813, 19.2338119002, 25.0463524709, 29.3718229816, 37.6100409447, 35.8947298854, 13.9427038167, 20.2988390265, 29.6773641371, 31.9110209761, 38.5532439532, 19.1488997281, 19.9281158996, 24.5916204474, 30.6975937394, 39.8329588338, 16.071110644, 24.7527151225, 32.3503618092, 37.332650551, 45.4405828945, 17.6319986283, 24.964777191, 24.2640949201, 34.7718166759, 37.1430784648, 18.1385210664, 23.0186156785, 30.8317753662, 35.5237516465, 37.7101699704, 13.668654211, 22.3640710385, 28.2924041402, 32.2187897206, 39.7382757559, 13.3069971355, 18.8612512195, 24.2412843677, 29.0166902632, 35.6948949398, 21.8627472096, 26.0927972681, 27.9554363863, 36.2551300861, 42.7344261141, 17.7399083381, 21.8350364849, 23.970998436, 29.0847355821, 39.7434142079, 11.3280462185, 22.1096931816, 25.731156431, 31.3480680122, 38.2880695717, 19.7734427391, 25.2340505623, 29.4650692613, 34.016337793, 39.4842497007, 18.0263843935, 24.1222412091, 33.7912164593, 31.5190672259, 42.1594426348, 18.6321295954, 25.8373692138, 33.2067650843, 39.2118721392, 38.5425202667, 15.6189433573, 18.3813002337, 26.7362005713, 29.7520066983, 33.8386141268, 14.2954900753, 21.8888575362, 34.7634255859, 36.6550895554, 42.7374615767, 18.5779737286, 23.2671614959, 32.1352681826, 36.5582063071, 41.344526313, 15.9045026683, 21.6786418167, 22.5542976518, 32.8611533201, 34.6477804281, 14.2962917155, 17.6511160154, 27.493741768, 32.6293618821, 34.8487613312, 15.9856491352, 19.4291180899, 22.9866068107, 33.572138553, 39.4527007301, 12.9783996202, 16.039845142, 25.5656458624, 30.2656397401, 37.8234451206, 18.0348586956, 21.7894172607, 25.539539217, 34.8870800713, 37.0497381227, 13.7665569953, 15.3372042678, 29.031712764, 36.2232424294, 41.389088496, 23.9093926265, 23.4016097939, 32.5022993206, 44.2368143493, 50.2476515212, 16.1771202166, 24.4283075251, 36.4818276343, 37.9615225989, 46.6396184266, 12.3106038912, 22.1608632252, 30.4740890289, 36.6387701863, 47.2135718405, 23.7521090996, 29.1071176962, 37.378026201, 39.2997196733, 48.5119199352, 16.7342482165, 19.9405124667, 32.4949213136, 34.6313825824, 44.6154715752, 13.7567821978, 19.4543947464, 33.4985631304, 36.971025909, 42.3967332642, 17.1116997022, 21.1477584542, 36.9712435205, 39.9442307239, 47.3207348393, 15.8532284677, 22.8075277, 35.09269263, 39.3256904945, 45.3533510015, 16.6393011109, 20.5897617204, 34.6254576978, 37.8026061809, 45.0104580842, 16.5267134308, 21.6687304172, 27.9470362026, 38.6357290375, 46.2211434503, 16.8891822501, 21.993528475, 31.9623264909, 39.2169270916, 42.6909489869, 15.1652197918, 21.6714127933, 31.8167411636, 40.5516317765, 46.1123835584, 13.284817119, 26.3654502694, 33.0205373355, 37.2048032491, 45.4630065794, 17.7928158611, 25.867915902, 32.391614397, 38.5982862225, 40.109818887, 15.5739329206, 24.7576854084, 29.9331546638, 36.3015808517, 42.2961668054, 19.1289228224, 24.5157482386, 26.7158840889, 38.5972013974, 45.044833904, 17.2630860403, 19.9160132486, 29.5096833409, 36.1992921424, 41.8353370405, 25.342137812, 27.8334238723, 35.1259923831, 47.2025870747, 48.6996888157, 10.0767341274, 25.5220099753, 23.3514811693, 32.9850542777, 38.3893962071, 18.649199321, 22.6754870957, 27.0380773442, 40.2922327962, 43.1097859763, 12.9569832879, 16.3383965688, 25.4558397479, 38.0585254191, 38.7359400597, 15.0369800453, 20.1527928049, 29.7112161784, 36.0745496239, 40.8362211606, 19.8629965066, 22.048112261, 30.2355044588, 40.7230413841, 45.7605291066, 13.1327405347, 18.9471370539, 25.8820678914, 35.8942078967, 43.7758349731, 20.3127370803, 21.2630571117, 27.3351814653, 34.9786303814, 48.4552254719, 15.7065455452, 22.1665111988, 27.6299744653, 32.3493109453, 37.0980958916, 15.9955551913, 20.9615466977, 23.0593149504, 29.6732562268, 40.0194944336, 11.9680800721, 13.7102182341, 23.8291144424, 29.3642212821, 38.2115456628, 12.5681391749, 18.3833793528, 24.6290189957, 31.4683766673, 38.6503274544, 15.5061636679, 17.5002447558, 24.7614128909, 32.9480510185, 36.8960058371, 20.3339266733, 22.9009052795, 30.4286001364, 34.4131902208, 38.2874849501, 14.4572499653, 18.3594472712, 27.0365426279, 33.8227932637, 37.7275957174, 15.3806331188, 24.0924243088, 29.9055973715, 34.4022978705, 34.3851568626, 16.0327690171, 22.9098092214, 31.9579692903, 32.8618924356, 45.3109341999, 12.1774573013, 23.079355808, 30.2079593127, 35.3724628459, 40.4560489198, 16.7717117177, 18.3397185105, 22.0330063514, 31.229723733, 36.6049602287, 13.6360769294, 22.4552041012, 28.0651009735, 36.2733183195, 35.0194230999, 17.7432322013, 26.2002298996, 29.350287087, 38.1087832329, 39.7599358968, 17.2279501593, 19.1620963295, 23.9250626673, 33.1185318786, 37.5763975976, 14.9021781849, 21.5987318339, 27.2507391452, 36.7900972862, 39.6807224568, 14.251778125, 22.4995294622, 27.7581312372, 29.6647105455, 37.1291492364, 12.8843584761, 24.099516675, 26.7532526424, 32.1502965412, 35.6721445267, 16.7135599501, 20.7949603911, 25.7243267394, 32.1454915883, 39.4472669133, 16.2824455578, 25.2389056242, 30.6441994802, 33.156334176, 41.8405747473, 14.5212496313, 20.8504996739, 26.63358806, 33.0640258056, 37.5476227201, 16.1305618056, 23.0091974471, 31.5278742505, 36.0428089639, 38.6675147675, 14.7036683143, 25.0324773284, 30.4220045509, 30.601759868, 36.4849708263, 18.369811284, 23.0410961966, 26.7564577303, 30.8158726843, 41.3456343761, 17.5033808961, 28.3634245767, 30.7413846451, 36.6750841385, 38.6122578974, 15.7837047065, 22.7304883852, 23.357527592, 35.1963583827, 42.0363147254] + }, + "_ddd_data": { + "obs_id": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200], + "group": [1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0], + "partition": [0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0], + "time": [0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1], + "outcome": [11.5190071583, 11.0668309674, 16.1968625766, 10.4436315868, 15.9593015248, 11.4484551777, 23.1671695492, 15.7263542986, 11.5321546753, 8.761747672, 11.4922379662, 13.199717822, 13.46650767, 16.7876752456, 12.8244741298, 8.9282176158, 15.6632068825, 12.6372615844, 13.590013548, 17.4324219277, 11.0073074889, 22.4546502976, 12.0848980587, 16.3955655823, 12.7702218611, 10.8366190685, 13.7549441387, 15.6894589213, 13.8647779785, 9.8492240111, 10.0509928699, 12.6430087, 22.4831938638, 15.4936443736, 12.1514558929, 15.4158910297, 16.3688067326, 15.7946543397, 15.7207406267, 10.6637633451, 12.2007488184, 22.5541966223, 21.1636934072, 13.905411838, 10.2049585806, 13.654912022, 22.2526117034, 20.7059975345, 9.0408295556, 17.0857748537, 16.4037749047, 13.5864875367, 17.3152284462, 13.1288214286, 7.9990707623, 13.3337771974, 17.1713251274, 14.0595392423, 14.1231384018, 10.8491444344, 12.2941786052, 14.4459442179, 10.8542562769, 13.8146220323, 8.7987779493, 14.036972167, 10.1077747449, 12.9158918995, 11.9956196416, 16.0374151861, 21.867911963, 15.4767874236, 15.2829697899, 11.7163977959, 10.3856678904, 15.6484871265, 21.4782039066, 14.4318687993, 11.9283659033, 21.8259817656, 10.5156677286, 21.7656347227, 10.8414965742, 14.2502366041, 15.7282362849, 16.4479519959, 14.7984175699, 13.5338839036, 9.7306486048, 11.1090345919, 12.848707012, 13.9772352987, 12.2442258511, 13.0576282921, 10.7707827235, 12.9980689086, 15.2584816646, 17.2488636888, 12.6193629505, 17.5499606936, 16.5168728298, 11.4732825359, 13.7036077788, 21.0286147708, 10.9038437584, 10.0490504509, 11.8015041434, 14.1900189986, 17.2977058996, 14.966126277, 13.2615592458, 15.5465639395, 14.9824038802, 12.6167160401, 14.8727554117, 22.969545014, 14.383846665, 8.1484443369, 11.4460032632, 14.0647732143, 10.8131950374, 15.3091835259, 19.3000701914, 14.0609666388, 12.0737516975, 12.0458035797, 12.1574125402, 12.4315653729, 13.6034502639, 15.3099782258, 15.9703933999, 8.7573297294, 17.3815754564, 15.204458937, 13.8240739637, 8.3373705978, 10.9306936564, 12.6355138173, 22.0437220076, 14.3480123037, 17.9595935489, 10.6816196756, 19.8867998851, 15.7736952724, 15.3124031588, 12.446041053, 13.1876152762, 15.2120554803, 15.8762940284, 11.022664494, 21.8337385085, 13.8625633836, 22.0973404852, 8.3743832608, 13.9953792322, 12.2602421677, 14.0389909129, 16.2350721416, 15.853527373, 21.9421126646, 16.4823694661, 12.9929436368, 14.753604502, 11.9665124752, 11.4290378188, 15.2410793463, 10.9656406391, 21.369268046, 10.58680772, 11.5836773438, 10.7151121905, 22.1634163195, 14.2632857649, 15.0458737762, 9.5154045838, 22.1891288117, 15.5510063316, 11.4997593311, 15.8093820424, 9.1746720429, 13.1454704522, 11.5315731876, 12.1647941947, 9.931236351, 22.746771687, 11.0744812656, 9.2279177648, 16.1527641067, 13.9885968452, 21.9265416653, 8.6129734464, 20.6933240956, 9.2316046749, 10.9728918746, 12.978572935, 16.170498071, 12.5653829614, 10.8861202167, 12.6071059949, 13.2754569687], + "stratum": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], + "psu": [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10, 8, 9, 10], + "fpc": [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180], + "weight": [2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.6666666667, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 2.7272727273, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] + }, + "_metadata": { + "description": "Golden values for survey-enabled estimator validation", + "r_version": "4.5.2", + "survey_version": "4.5", + "seed": 42, + "n_units": 150, + "n_periods": 5, + "n_obs_staggered": 750, + "n_obs_ddd": 200 + } +} diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index c1fa24ce..2943eff2 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -746,8 +746,84 @@ comparisons match to machine precision (differences < 1e-10). ``subset()`` drops empty strata (df=199). This is a documented deviation (see REGISTRY.md); the diff-diff approach is conservative per Lumley (2004). -Reproducing Survey Validation -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Survey Estimator Validation +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Four additional estimators are validated against R's ``survey::svyglm()`` using +synthetic staggered-adoption and DDD datasets. Each estimator reduces to a WLS +regression under survey weights, so the R comparison fits the equivalent +``svyglm()`` model and compares coefficients and standard errors. + +**Data:** 150-unit staggered panel (5 periods, 4 strata, 10 PSUs, FPC) with +cohorts at *t* = 3 and *t* = 4; 200-observation DDD cross-section (4 strata, +10 PSUs, FPC). Both generated with seed 42. + +.. list-table:: + :header-rows: 1 + :widths: 8 18 30 15 15 14 + + * - Test + - Estimator + - R Comparison + - Coef Gap + - SE Gap + - Tolerance + * - S1 + - ``ImputationDiD`` + - ``svyglm()`` on control-only (Omega_0) FE regression; covariate coefficients + - < 1e-10 + - 0.00% + - 1.5% + * - S2 + - ``StackedDiD`` + - ``svyglm()`` on stacked dataset with Q-weight x survey weight composition + - < 1e-10 + - 0.77% + - 1.5% + * - S3 + - ``SunAbraham`` + - ``svyglm()`` with cohort x period interactions; IW-aggregated ATT + - < 1e-11 + - 0.00% + - 1.5% + * - S4 + - ``TripleDifference`` + - ``svyglm()`` three-way interaction (``group:partition:time``) + - < 1e-10 + - 0.36% + - 1.5% + +**Key details:** + +- **S1** validates the WLS building block that ``ImputationDiD`` uses internally + (control-only regression with absorbed unit + time FE and time-varying + covariates). A companion smoke test confirms ``ImputationDiD.fit()`` produces + finite ATT/SE under survey weights. +- **S2** replicates the full stacking pipeline in R: sub-experiment construction, + sample-share Q-weight computation, Q x survey weight composition with + normalization, then ``svyglm()`` on the stacked data with strata/PSU structure. + The 0.77% SE gap arises because R omits FPC on the stacked data while Python + re-resolves the full survey design. +- **S3** compares both individual cohort x relative-time effects and the + IW-aggregated overall ATT (with survey-weighted cohort masses and delta-method + SE via the vcov submatrix). +- **S4** exploits the algebraic equivalence between the pairwise DDD + decomposition (``estimation_method="reg"``, no covariates) and the three-way + interaction coefficient from a single OLS regression. + +Reproducing Survey Estimator Validation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: bash + + # Generate golden values + Rscript benchmarks/R/benchmark_survey_estimators.R + + # Run validation tests + pytest tests/test_survey_estimator_validation.py -v + +Reproducing Survey Real-Data Validation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: bash diff --git a/tests/test_survey_estimator_validation.py b/tests/test_survey_estimator_validation.py new file mode 100644 index 00000000..fed08046 --- /dev/null +++ b/tests/test_survey_estimator_validation.py @@ -0,0 +1,325 @@ +"""Cross-validation tests: survey-enabled estimators vs R's survey::svyglm(). + +Validates 4 estimators against R golden values: + S1: ImputationDiD — control-only WLS regression (covariate coefficients) + S2: StackedDiD — stacked WLS with Q-weight x survey weight composition + S3: SunAbraham — interaction-weighted ATT + S4: TripleDifference — three-way interaction DDD + +JSON fixtures generated by: benchmarks/R/benchmark_survey_estimators.R +""" + +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from diff_diff import SurveyDesign +from diff_diff.imputation import ImputationDiD +from diff_diff.linalg import LinearRegression +from diff_diff.stacked_did import StackedDiD +from diff_diff.sun_abraham import SunAbraham +from diff_diff.triple_diff import TripleDifference + +GOLDEN_FILE = ( + Path(__file__).parent.parent + / "benchmarks" + / "data" + / "synthetic" + / "survey_estimator_validation_golden.json" +) + +# Tolerances — observed gaps: S1 0.00%, S2 0.77%, S3 0.00%, S4 0.36% +COEF_RTOL = 1e-8 +COEF_ATOL = 1e-8 +SE_RTOL = 0.015 # 1.5% — covers S2 (0.77%) and S4 (0.36%) with headroom +SE_ATOL = 1e-6 + + +def _assert_close(py_val, r_val, rtol, atol, label): + """Assert values are close with informative error message.""" + if np.isnan(r_val): + return + diff = abs(py_val - r_val) + threshold = max(atol, rtol * abs(r_val)) + assert diff < threshold, ( + f"{label}: Python={py_val:.6f}, R={r_val:.6f}, " + f"diff={diff:.6f}, threshold={threshold:.6f} " + f"(rtol={rtol}, atol={atol})" + ) + + +def _load_embedded_data(results, key): + """Load embedded dataset from golden JSON.""" + data = results[key] + return pd.DataFrame({col: vals for col, vals in data.items()}) + + +def _parse_cohort_key(key_str): + """Parse 'g3_e0' -> (3, 0), 'g3_em2' -> (3, -2).""" + g_part, e_part = key_str.split("_") + g = int(g_part[1:]) # 'g3' -> 3 + e_str = e_part[1:] # 'e0' -> '0', 'em2' -> 'm2' + e = -int(e_str[1:]) if e_str.startswith("m") else int(e_str) + return g, e + + +@pytest.fixture(scope="module") +def golden(): + """Load pre-computed R golden values.""" + if not GOLDEN_FILE.exists(): + pytest.skip( + f"Golden values not found: {GOLDEN_FILE}. " + "Run: Rscript benchmarks/R/benchmark_survey_estimators.R" + ) + with open(GOLDEN_FILE) as f: + return json.load(f) + + +class TestSurveyEstimatorValidation: + """Validate survey-enabled estimators against R's survey::svyglm().""" + + def _load_staggered_data(self, golden): + """Load staggered panel from _staggered_data key.""" + df = _load_embedded_data(golden, "_staggered_data") + # Handle Inf serialization: R's toJSON writes Inf as "Inf" (string) + df["first_treat"] = pd.to_numeric(df["first_treat"], errors="coerce") + for col in ["unit", "period"]: + df[col] = df[col].astype(int) + for col in ["outcome", "weight", "fpc", "x1", "x2"]: + df[col] = pd.to_numeric(df[col], errors="coerce") + df["stratum"] = df["stratum"].astype(int) + df["psu"] = df["psu"].astype(int) + return df + + def _load_ddd_data(self, golden): + """Load DDD cross-section from _ddd_data key.""" + df = _load_embedded_data(golden, "_ddd_data") + for col in ["group", "partition", "time"]: + df[col] = df[col].astype(int) + for col in ["outcome", "weight", "fpc"]: + df[col] = pd.to_numeric(df[col], errors="coerce") + df["stratum"] = df["stratum"].astype(int) + df["psu"] = df["psu"].astype(int) + return df + + # ------------------------------------------------------------------ + # S1: ImputationDiD — control-only WLS regression + # ------------------------------------------------------------------ + + def test_s1_imputation_control_regression(self, golden): + """Validate control-only WLS covariate coefficients against R's svyglm.""" + r = golden["s1_imputation_did"] + data = self._load_staggered_data(golden) + + # Build Omega_0 (untreated observations) + # Matches imputation.py line 341: omega_0_mask = ~df["_treated"] + omega_0 = data[ + (np.isinf(data["first_treat"])) + | (data["period"] < data["first_treat"]) + ].copy() + omega_0 = omega_0.reset_index(drop=True) + + sd = SurveyDesign( + weights="weight", strata="stratum", psu="psu", fpc="fpc" + ) + resolved = sd.resolve(omega_0) + survey_weights = omega_0["weight"].values.astype(float) + + # Build design matrix: unit + time dummies + covariates + # Algebraically equivalent to within-transformation for covariate coefs + unit_dummies = pd.get_dummies( + omega_0["unit"].astype(str), prefix="u", drop_first=True + ) + time_dummies = pd.get_dummies( + omega_0["period"].astype(str), prefix="t", drop_first=True + ) + X = np.column_stack( + [ + unit_dummies.values.astype(float), + time_dummies.values.astype(float), + omega_0["x1"].values, + omega_0["x2"].values, + ] + ) + y = omega_0["outcome"].values + + reg = LinearRegression( + include_intercept=True, + survey_design=resolved, + weights=survey_weights, + weight_type="pweight", + ) + reg.fit(X, y) + + # x1 and x2 are last two columns + _assert_close( + reg.coefficients_[-2], + r["coef_x1"], + COEF_RTOL, + COEF_ATOL, + "S1 x1 coef", + ) + _assert_close( + reg.coefficients_[-1], + r["coef_x2"], + COEF_RTOL, + COEF_ATOL, + "S1 x2 coef", + ) + _assert_close( + np.sqrt(reg.vcov_[-2, -2]), + r["se_x1"], + SE_RTOL, + SE_ATOL, + "S1 x1 SE", + ) + _assert_close( + np.sqrt(reg.vcov_[-1, -1]), + r["se_x2"], + SE_RTOL, + SE_ATOL, + "S1 x2 SE", + ) + + def test_s1_imputation_end_to_end(self, golden): + """Smoke test: ImputationDiD.fit() with survey weights produces finite ATT.""" + data = self._load_staggered_data(golden) + sd = SurveyDesign( + weights="weight", strata="stratum", psu="psu", fpc="fpc" + ) + est = ImputationDiD() + result = est.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["x1", "x2"], + survey_design=sd, + ) + assert np.isfinite(result.overall_att), "ImputationDiD ATT should be finite" + assert np.isfinite(result.overall_se), "ImputationDiD SE should be finite" + assert result.overall_se > 0, "ImputationDiD SE should be positive" + + # ------------------------------------------------------------------ + # S2: StackedDiD — stacked WLS with Q-weight x survey weight + # ------------------------------------------------------------------ + + def test_s2_stacked_did(self, golden): + """Validate StackedDiD ATT against R's svyglm on stacked data.""" + r = golden["s2_stacked_did"] + data = self._load_staggered_data(golden) + + sd = SurveyDesign( + weights="weight", strata="stratum", psu="psu", fpc="fpc" + ) + est = StackedDiD( + kappa_pre=1, + kappa_post=1, + weighting="sample_share", + clean_control="never_treated", + ) + result = est.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="event_study", + survey_design=sd, + ) + + _assert_close( + result.overall_att, r["att"], COEF_RTOL, COEF_ATOL, "S2 ATT" + ) + _assert_close( + result.overall_se, r["se"], SE_RTOL, SE_ATOL, "S2 SE" + ) + + # Compare individual event-study coefficients + if "event_study" in r: + for e_str, r_eff in r["event_study"].items(): + e = int(e_str) + py_eff = result.event_study_effects.get(e) + if py_eff is not None: + _assert_close( + py_eff["effect"], + r_eff["coef"], + COEF_RTOL, + COEF_ATOL, + f"S2 delta_{e}", + ) + + # ------------------------------------------------------------------ + # S3: SunAbraham — interaction-weighted ATT + # ------------------------------------------------------------------ + + def test_s3_sun_abraham(self, golden): + """Validate SunAbraham IW-aggregated ATT against R's svyglm.""" + r = golden["s3_sun_abraham"] + data = self._load_staggered_data(golden) + + sd = SurveyDesign( + weights="weight", strata="stratum", psu="psu", fpc="fpc" + ) + est = SunAbraham() + result = est.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + survey_design=sd, + ) + + # Compare overall IW-aggregated ATT + _assert_close( + result.overall_att, r["att"], COEF_RTOL, COEF_ATOL, "S3 ATT" + ) + _assert_close( + result.overall_se, r["se"], SE_RTOL, SE_ATOL, "S3 SE" + ) + + # Compare individual cohort x rel-time effects (post-treatment only) + if "cohort_effects" in r: + for key_str, r_delta in r["cohort_effects"].items(): + g, e = _parse_cohort_key(key_str) + py_entry = result.cohort_effects.get((g, e)) + if py_entry is not None: + _assert_close( + py_entry["effect"], + r_delta["coef"], + COEF_RTOL, + COEF_ATOL, + f"S3 delta_{g}_{e}", + ) + + # ------------------------------------------------------------------ + # S4: TripleDifference — three-way interaction DDD + # ------------------------------------------------------------------ + + def test_s4_triple_diff(self, golden): + """Validate TripleDifference DDD against R's svyglm three-way interaction.""" + r = golden["s4_triple_diff"] + data = self._load_ddd_data(golden) + + sd = SurveyDesign( + weights="weight", strata="stratum", psu="psu", fpc="fpc" + ) + est = TripleDifference(estimation_method="reg") + result = est.fit( + data, + outcome="outcome", + group="group", + partition="partition", + time="time", + survey_design=sd, + ) + + _assert_close( + result.att, r["att"], COEF_RTOL, COEF_ATOL, "S4 DDD coef" + ) + _assert_close(result.se, r["se"], SE_RTOL, SE_ATOL, "S4 DDD SE")