Batch Cohort

#!/usr/bin/env Rscript
# ==============================================================================
# KNHANES Survey-Weighted Logistic Regression — Batch Template
# Combo: {{COMBO_ID}} | {{EXPOSURE_LABEL}} → {{OUTCOME_LABEL}}
#
# Slot variables (replaced by batch_template_generator.R):
#   {{EXPOSURE_VAR}}, {{EXPOSURE_LABEL}}, {{EXPOSURE_CODING}}
#   {{OUTCOME_VAR}}, {{OUTCOME_LABEL}}, {{OUTCOME_CODING}}
#   {{RESULTS_DIR}}, {{COMBO_ID}}, {{COVARIATES_REMOVE}}
#
# Data: KNHANES single-cycle CSV (e.g., HN18.csv)
# ==============================================================================

suppressPackageStartupMessages({
  library(survey)
  library(dplyr)
  library(tidyr)
})

# --- CONFIG (edit these paths for your setup) ---
DATA_PATH     <- "{{DATA_PATH}}"       # e.g., /path/to/HN18.csv
RESULTS_DIR   <- "{{RESULTS_DIR}}"
COMBO_ID      <- "{{COMBO_ID}}"
AGE_MIN       <- 20
AGE_MAX       <- Inf

# --- Standard Covariates ---
# Full set; exposure-specific removals applied below
ALL_COVARIATES <- c("age", "sex_binary", "edu_binary", "income_binary",
                     "smoking_3cat", "alcohol_3cat", "obesity_binary", "cvd_binary")
REMOVE_COVARS  <- c({{COVARIATES_REMOVE}})
COVARIATES     <- setdiff(ALL_COVARIATES, REMOVE_COVARS)

# --- Load & Clean ---
raw <- read.csv(DATA_PATH)

df <- raw %>%
  filter(age >= AGE_MIN, age <= AGE_MAX) %>%
  mutate(
    # --- Demographics ---
    sex_binary     = ifelse(sex == 1, 0, 1),  # 0=male, 1=female
    edu_binary     = ifelse(edu %in% 1:3, 0, 1),  # 0=non-college, 1=college
    income_binary  = ifelse(incm %in% 1:3, 0, 1),  # 0=bottom80, 1=top20

    # --- Standard Covariates ---
    smoking_3cat = case_when(
      BS3_1 %in% c(1, 2) ~ "current",
      BS3_1 == 3          ~ "former",
      BS3_1 == 8          ~ "never",
      TRUE                ~ NA_character_
    ),
    alcohol_3cat = case_when(
      BD1_11 %in% 2:6 ~ "frequent",
      BD1_11 == 1      ~ "occasional",
      BD1_11 == 8      ~ "never",
      TRUE             ~ NA_character_
    ),
    obesity_binary = ifelse(HE_obe >= 4, 1, 0),
    cvd_binary     = ifelse(DI4_dg == 1 | DI5_dg == 1 | DI6_dg == 1, 1, 0),

    # --- Exposure (SLOT) ---
    exposure = {{EXPOSURE_CODING}},

    # --- Outcome (SLOT) ---
    outcome = {{OUTCOME_CODING}}
  ) %>%
  filter(!is.na(exposure), !is.na(outcome))

cat(sprintf("[%s] N after filtering: %d (exposure=1: %d, outcome=1: %d)\n",
            COMBO_ID, nrow(df), sum(df$exposure == 1), sum(df$outcome == 1)))

# --- EPV Check ---
n_events <- min(sum(df$outcome == 1), sum(df$outcome == 0))
n_params <- length(COVARIATES) + 1  # +1 for exposure
epv <- n_events / n_params
if (epv < 10) {
  warning(sprintf("[%s] EPV = %.1f (< 10). Results may be unreliable.", COMBO_ID, epv))
}

# --- Survey Design ---
design <- svydesign(
  id      = ~psu,
  strata  = ~kstrata,
  weights = ~wt_itvex,
  nest    = TRUE,
  data    = df
)

# --- Table 1: Weighted Demographics by Exposure ---
table1_vars <- c("age", "sex_binary", "edu_binary", "income_binary",
                 "smoking_3cat", "alcohol_3cat", "obesity_binary", "cvd_binary", "outcome")

table1_rows <- lapply(table1_vars, function(v) {
  tryCatch({
    if (is.numeric(df[[v]])) {
      means <- svyby(as.formula(paste0("~", v)), ~exposure, design, svymean, na.rm = TRUE)
      data.frame(variable = v,
                 unexposed_mean = means[means$exposure == 0, v],
                 exposed_mean   = means[means$exposure == 1, v],
                 stringsAsFactors = FALSE)
    } else {
      props <- svyby(as.formula(paste0("~factor(", v, ")")), ~exposure, design, svymean, na.rm = TRUE)
      data.frame(variable = paste0(v, "_", names(props)[-1]),
                 unexposed_mean = as.numeric(props[props$exposure == 0, -1]),
                 exposed_mean   = as.numeric(props[props$exposure == 1, -1]),
                 stringsAsFactors = FALSE)
    }
  }, error = function(e) NULL)
})
table1 <- do.call(rbind, Filter(Negate(is.null), table1_rows))
write.csv(table1, file.path(RESULTS_DIR, "table1.csv"), row.names = FALSE)

# --- Helper: extract OR + CI + p from svyglm ---
extract_or <- function(model) {
  or_val  <- unname(exp(coef(model)["exposure"]))
  ci_vals <- unname(exp(confint(model))["exposure", ])
  p_val   <- unname(summary(model)$coefficients["exposure", "Pr(>|t|)"])
  c(OR = or_val, lower = ci_vals[1], upper = ci_vals[2], p = p_val)
}

# --- Main Analysis: Sequential Models ---
results <- list()

# Model 1: Unadjusted
m1 <- svyglm(outcome ~ exposure, design = design, family = quasibinomial())
results[["Model 1 (unadjusted)"]] <- extract_or(m1)

# Model 2: Age + Sex
m2 <- svyglm(outcome ~ exposure + age + sex_binary, design = design, family = quasibinomial())
results[["Model 2 (age+sex)"]] <- extract_or(m2)

# Model 3: Fully adjusted
formula_full <- as.formula(paste("outcome ~ exposure +", paste(COVARIATES, collapse = " + ")))
m3 <- svyglm(formula_full, design = design, family = quasibinomial())
results[["Model 3 (fully adjusted)"]] <- extract_or(m3)

# Save main results
main_df <- data.frame(
  combo_id = COMBO_ID,
  exposure = "{{EXPOSURE_LABEL}}",
  outcome  = "{{OUTCOME_LABEL}}",
  model    = names(results),
  N        = nrow(df),
  events   = sum(df$outcome == 1),
  OR       = sapply(results, `[`, "OR"),
  lower    = sapply(results, `[`, "lower"),
  upper    = sapply(results, `[`, "upper"),
  p_value  = sapply(results, `[`, "p"),
  row.names = NULL
)
write.csv(main_df, file.path(RESULTS_DIR, "main_results.csv"), row.names = FALSE)

# --- Subgroup Analyses (Model 3 within strata) ---
subgroup_vars <- list(
  sex     = list(var = "sex_binary", levels = c(0, 1), labels = c("Male", "Female")),
  age_grp = list(var = "age_group",  levels = c("20-39", "40-59", "60+"), labels = NULL),
  edu     = list(var = "edu_binary", levels = c(0, 1), labels = c("Non-college", "College")),
  income  = list(var = "income_binary", levels = c(0, 1), labels = c("Lower", "Higher"))
)

# Create age groups
df$age_group <- cut(df$age, breaks = c(20, 40, 60, Inf), right = FALSE,
                    labels = c("20-39", "40-59", "60+"))
design <- update(design, age_group = df$age_group)

subgroup_results <- list()
for (sg_name in names(subgroup_vars)) {
  sg <- subgroup_vars[[sg_name]]
  for (lev in sg$levels) {
    tryCatch({
      sub_design <- subset(design, get(sg$var) == lev)
      n_sub <- nrow(sub_design$variables)
      n_events_sub <- sum(sub_design$variables$outcome == 1)

      if (n_sub >= 30 && n_events_sub >= 5) {
        # Reduced model for subgroups (fewer covariates to avoid convergence issues)
        sg_covars <- setdiff(COVARIATES, sg$var)
        sg_formula <- as.formula(paste("outcome ~ exposure +", paste(sg_covars, collapse = " + ")))
        m_sg <- svyglm(sg_formula, design = sub_design, family = quasibinomial())
        r_sg <- extract_or(m_sg)

        subgroup_results[[length(subgroup_results) + 1]] <- data.frame(
          combo_id  = COMBO_ID,
          exposure  = "{{EXPOSURE_LABEL}}",
          outcome   = "{{OUTCOME_LABEL}}",
          subgroup  = sg_name,
          level     = as.character(lev),
          N         = n_sub,
          events    = n_events_sub,
          OR        = r_sg["OR"],
          lower     = r_sg["lower"],
          upper     = r_sg["upper"],
          p_value   = r_sg["p"],
          row.names = NULL
        )
      }
    }, error = function(e) {
      message(sprintf("[%s] Subgroup %s=%s failed: %s", COMBO_ID, sg_name, lev, e$message))
    })
  }
}

if (length(subgroup_results) > 0) {
  subgroup_df <- do.call(rbind, subgroup_results)
  write.csv(subgroup_df, file.path(RESULTS_DIR, "subgroup_results.csv"), row.names = FALSE)
}

cat(sprintf("[%s] Complete. Fully adjusted wOR: %.2f (%.2f–%.2f), p=%.3f\n",
            COMBO_ID,
            main_df$OR[3], main_df$lower[3], main_df$upper[3], main_df$p_value[3]))

#!/usr/bin/env Rscript
# ==============================================================================
# Batch Cohort Analysis Generator
# Generates N analysis scripts from a base template × variable combinations
#
# Usage:
#   Rscript batch_template_generator.R \
#     --template base_template.R \
#     --combinations combinations.csv \
#     --outdir batch_20260417/ \
#     --mode code_only|execute|full
#
# Combinations CSV columns:
#   id, exposure_var, exposure_label, exposure_coding,
#   outcome_var, outcome_label, outcome_coding,
#   covariates_remove (comma-separated vars to drop from default set)
# ==============================================================================

suppressPackageStartupMessages({
  library(optparse)
  library(glue)
})

# --- CLI Arguments ---
option_list <- list(
  make_option("--template", type = "character", help = "Path to base R template"),
  make_option("--combinations", type = "character", help = "CSV with variable combinations"),
  make_option("--outdir", type = "character", default = "batch_output", help = "Output directory"),
  make_option("--mode", type = "character", default = "code_only",
              help = "code_only | execute | full"),
  make_option("--parallel", type = "integer", default = 1, help = "Number of parallel workers")
)
opt <- parse_args(OptionParser(option_list = option_list))

# --- Read Inputs ---
template_code <- readLines(opt$template)
combos <- read.csv(opt$combinations, stringsAsFactors = FALSE)
cat(glue("Loaded {nrow(combos)} combinations from {opt$combinations}\n\n"))

# --- Create Output Structure ---
dir.create(file.path(opt$outdir, "scripts"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(opt$outdir, "template"), showWarnings = FALSE)
dir.create(file.path(opt$outdir, "logs"), showWarnings = FALSE)

# Freeze template
writeLines(template_code, file.path(opt$outdir, "template", "base_template.R"))

# --- Generate Scripts ---
script_paths <- character(nrow(combos))

for (i in seq_len(nrow(combos))) {
  combo <- combos[i, ]
  script_name <- sprintf("%02d_%s_%s.R", i, combo$exposure_var, combo$outcome_var)
  script_path <- file.path(opt$outdir, "scripts", script_name)

  # Results directory for this combination
  results_dir <- file.path(opt$outdir, "results",
                           sprintf("%02d_%s_%s", i, combo$exposure_var, combo$outcome_var))
  dir.create(results_dir, recursive = TRUE, showWarnings = FALSE)

  # Replace slot variables in template
  code <- template_code
  code <- gsub("\\{\\{EXPOSURE_VAR\\}\\}", combo$exposure_var, code)
  code <- gsub("\\{\\{EXPOSURE_LABEL\\}\\}", combo$exposure_label, code)
  code <- gsub("\\{\\{EXPOSURE_CODING\\}\\}", combo$exposure_coding, code)
  code <- gsub("\\{\\{OUTCOME_VAR\\}\\}", combo$outcome_var, code)
  code <- gsub("\\{\\{OUTCOME_LABEL\\}\\}", combo$outcome_label, code)
  code <- gsub("\\{\\{OUTCOME_CODING\\}\\}", combo$outcome_coding, code)
  code <- gsub("\\{\\{RESULTS_DIR\\}\\}", results_dir, code)
  code <- gsub("\\{\\{COMBO_ID\\}\\}", sprintf("%02d", i), code)

  # Handle covariate removal (self-adjustment prevention)
  if (!is.na(combo$covariates_remove) && nchar(combo$covariates_remove) > 0) {
    remove_vars <- trimws(strsplit(combo$covariates_remove, ",")[[1]])
    remove_str <- paste0('"', remove_vars, '"', collapse = ", ")
    code <- gsub("\\{\\{COVARIATES_REMOVE\\}\\}", remove_str, code)
  } else {
    code <- gsub("\\{\\{COVARIATES_REMOVE\\}\\}", "", code)
  }

  writeLines(code, script_path)
  script_paths[i] <- script_path
  cat(glue("  Generated: {script_name}\n"))
}

# --- Generate Master Runner ---
runner_path <- file.path(opt$outdir, "scripts", "run_all.R")
runner_code <- c(
  "#!/usr/bin/env Rscript",
  "# Master runner — executes all batch scripts sequentially",
  "# Auto-generated by batch_template_generator.R",
  "",
  'cat("=== Batch Cohort Analysis ===\\n")',
  glue('cat("Combinations: {nrow(combos)}\\n")'),
  glue('cat("Started: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "\\n\\n")'),
  "",
  "results_log <- data.frame(",
  "  id = character(), exposure = character(), outcome = character(),",
  "  status = character(), time_sec = numeric(), error_msg = character(),",
  "  stringsAsFactors = FALSE",
  ")",
  "",
  "scripts <- c(",
  paste0('  "', script_paths, '"', collapse = ",\n"),
  ")",
  "",
  "for (i in seq_along(scripts)) {",
  "  cat(sprintf('\\n--- Running %d/%d: %s ---\\n', i, length(scripts), basename(scripts[i])))",
  "  t0 <- proc.time()['elapsed']",
  "  status <- tryCatch({",
  "    source(scripts[i], local = new.env())",
  '    "SUCCESS"',
  "  }, error = function(e) {",
  '    paste0("FAILED: ", conditionMessage(e))',
  "  })",
  "  elapsed <- round(proc.time()['elapsed'] - t0, 1)",
  "  is_ok <- startsWith(status, 'SUCCESS')",
  "  results_log <- rbind(results_log, data.frame(",
  paste0("    id = sprintf('%02d', i),"),
  "    exposure = basename(scripts[i]),",
  "    outcome = '',",
  '    status = ifelse(is_ok, "SUCCESS", "FAILED"),',
  "    time_sec = elapsed,",
  '    error_msg = ifelse(is_ok, "", status),',
  "    stringsAsFactors = FALSE",
  "  ))",
  "  cat(sprintf('  %s (%.1fs)\\n', status, elapsed))",
  "}",
  "",
  "# Save execution log",
  glue('write.csv(results_log, "{file.path(opt$outdir, "logs", "batch_execution.csv")}", row.names = FALSE)'),
  "",
  'cat("\\n=== Batch Complete ===\\n")',
  'cat(sprintf("Success: %d / %d\\n", sum(results_log$status == "SUCCESS"), nrow(results_log)))',
  'cat(sprintf("Failed:  %d / %d\\n", sum(results_log$status == "FAILED"), nrow(results_log)))'
)
writeLines(runner_code, runner_path)
cat(glue("\nMaster runner: {runner_path}\n"))

# --- Generate Summary Aggregator ---
aggregator_path <- file.path(opt$outdir, "scripts", "aggregate_results.R")
aggregator_code <- c(
  "#!/usr/bin/env Rscript",
  "# Aggregate results from all batch runs into summary matrix",
  "",
  glue('results_base <- "{file.path(opt$outdir, "results")}"'),
  "dirs <- list.dirs(results_base, recursive = FALSE)",
  "",
  "all_main <- list()",
  "all_subgroup <- list()",
  "",
  "for (d in dirs) {",
  '  main_file <- file.path(d, "main_results.csv")',
  "  if (file.exists(main_file)) {",
  "    df <- read.csv(main_file, stringsAsFactors = FALSE)",
  "    df$combo_id <- basename(d)",
  "    all_main[[length(all_main) + 1]] <- df",
  "  }",
  '  sub_file <- file.path(d, "subgroup_results.csv")',
  "  if (file.exists(sub_file)) {",
  "    df <- read.csv(sub_file, stringsAsFactors = FALSE)",
  "    df$combo_id <- basename(d)",
  "    all_subgroup[[length(all_subgroup) + 1]] <- df",
  "  }",
  "}",
  "",
  "if (length(all_main) > 0) {",
  "  summary_matrix <- do.call(rbind, all_main)",
  glue('  write.csv(summary_matrix, "{file.path(opt$outdir, "summary", "summary_matrix.csv")}", row.names = FALSE)'),
  '  cat(sprintf("Summary matrix: %d rows from %d combinations\\n", nrow(summary_matrix), length(all_main)))',
  "}",
  "",
  "if (length(all_subgroup) > 0) {",
  "  subgroup_matrix <- do.call(rbind, all_subgroup)",
  glue('  write.csv(subgroup_matrix, "{file.path(opt$outdir, "summary", "subgroup_matrix.csv")}", row.names = FALSE)'),
  '  cat(sprintf("Subgroup matrix: %d rows\\n", nrow(subgroup_matrix)))',
  "}"
)
dir.create(file.path(opt$outdir, "summary"), showWarnings = FALSE)
writeLines(aggregator_code, aggregator_path)
cat(glue("Aggregator: {aggregator_path}\n"))

# --- Execute if requested ---
if (opt$mode %in% c("execute", "full")) {
  cat("\n=== Executing batch ===\n")
  source(runner_path)

  cat("\n=== Aggregating results ===\n")
  source(aggregator_path)
}

# --- Generate README ---
readme_path <- file.path(opt$outdir, "README.md")
readme <- c(
  glue("# Batch Cohort Analysis — {Sys.Date()}"),
  "",
  glue("- **Template**: `{opt$template}`"),
  glue("- **Combinations**: {nrow(combos)}"),
  glue("- **Mode**: {opt$mode}"),
  "",
  "## Combination Matrix",
  "",
  "| # | Exposure | Outcome | Covariates Removed |",
  "|---|----------|---------|-------------------|",
  apply(combos, 1, function(row) {
    sprintf("| %s | %s | %s | %s |",
            row["id"], row["exposure_label"], row["outcome_label"],
            ifelse(is.na(row["covariates_remove"]), "—", row["covariates_remove"]))
  }),
  "",
  "## Notes",
  "",
  "- This is an **exploratory** analysis for hypothesis generation, not confirmation.",
  "- Multiple comparisons: interpret p-values with caution (Bonferroni threshold shown in summary).",
  "- Survey weights applied to all analyses.",
  glue("- Data file SHA256: [compute with `shasum -a 256 <datafile>`]"),
  ""
)
writeLines(readme, readme_path)

cat(glue("\n✅ Batch generation complete: {opt$outdir}/\n"))
cat(glue("   {nrow(combos)} scripts generated\n"))

schema_version: 2
name: batch-cohort
layer: B
owner_domain: batch_analysis
maturity: official

when_to_use: "Generate N analysis scripts from one validated methodology template across many exposure/outcome combinations."
when_NOT_to_use: "A single analysis (use analyze-stats); cross-country comparison (use cross-national)."

inputs:
  - "methodology template"
  - "exposure/outcome combination matrix"
outputs:
  - "N analysis scripts"
  - "summary matrix"
side_effects:
  - writes_project_artifacts
downstream_consumers:
  - analyze-stats
  - self-review
forbidden_actions:
  - alter_methodology_per_combination_without_disclosure
  - fabricate_results_matrix

# v2.1 quality card
purpose: "Scale one validated method across many variable combinations, swapping only exposure/outcome, never the method."
safety_boundaries:
  - "The methodology is held constant across all generated scripts; deviations are disclosed."
  - "Generated scripts run on real data; no results are pre-filled."
known_limitations:
  - "Inherits the source template's assumptions; a flawed template propagates."
  - "No standalone demo; outputs are code to be executed and reviewed."
validation_commands:
  - "execute each generated script and reconcile the summary matrix"
  - "/self-review"
evidence_surface: manual_workflow