Data → Validation → Preparation → Analysis → Diagnostics → Reporting

# Standard function signature pattern
itc_function(
  data,                    # Primary data input
  outcome_var,             # Outcome variable name
  treatment_var,           # Treatment variable name
  covariates = NULL,       # Optional covariates
  method = "default",      # Method specification
  alpha = 0.05,           # Significance level
  seed = NULL,            # For reproducibility
  verbose = TRUE,         # Progress messages
  ...                      # Additional method-specific args
)

# Standard return structure
list(
  results = tibble(...),   # Main results as tibble
  diagnostics = list(...), # Model diagnostics
  model = fitted_model,    # Raw model object
  data_summary = list(...),# Data summary
  call = match.call(),     # Original call
  parameters = list(...)   # Analysis parameters
)

# Recommended project structure
project/
├── R/
│   ├── 01_data_prep.R
│   ├── 02_analysis.R
│   ├── 03_sensitivity.R
│   └── 04_reporting.R
├── data/
│   ├── raw/
│   └── processed/
├── output/
│   ├── figures/
│   └── tables/
├── renv.lock          # Package versions
└── _targets.R         # Pipeline definition (optional)

# Load packages with explicit namespacing preference
library(tidyverse)
library(meta)       # Pairwise MA
library(netmeta)    # NMA
library(maicplus)   # MAIC
library(stc)        # STC
library(multinma)   # ML-NMR

# Set global options
options(
  dplyr.summarise.inform = FALSE,
  mc.cores = parallel::detectCores() - 1
)

# Set seed for reproducibility
set.seed(12345)

# Validate IPD structure
validate_ipd <- function(data, outcome_var, treatment_var, covariates = NULL) {
  errors <- character()
  warnings <- character()


  # Check required columns exist
  required_cols <- c(outcome_var, treatment_var)
  if (!is.null(covariates)) required_cols <- c(required_cols, covariates)

  missing_cols <- setdiff(required_cols, names(data))
  if (length(missing_cols) > 0) {
    errors <- c(errors, paste("Missing columns:", paste(missing_cols, collapse = ", ")))
  }

  # Check outcome type
  if (outcome_var %in% names(data)) {
    outcome_vals <- unique(data[[outcome_var]])
    if (all(outcome_vals %in% c(0, 1, NA))) {
      message("Detected binary outcome")
    } else if (is.numeric(data[[outcome_var]])) {
      message("Detected continuous outcome")
    }
  }

  # Check treatment levels
  if (treatment_var %in% names(data)) {
    n_trt <- length(unique(data[[treatment_var]]))
    if (n_trt < 2) {
      errors <- c(errors, "Treatment variable must have at least 2 levels")
    }
    message(sprintf("Found %d treatment levels", n_trt))
  }

  # Check for missing values
  if (any(is.na(data[required_cols]))) {
    n_missing <- sum(!complete.cases(data[required_cols]))
    warnings <- c(warnings, sprintf("%d observations with missing values", n_missing))
  }

  list(
    valid = length(errors) == 0,
    errors = errors,
    warnings = warnings,
    n_obs = nrow(data),
    n_complete = sum(complete.cases(data[required_cols]))
  )
}

# Create preparation recipe
create_itc_recipe <- function(data, outcome_var, treatment_var, covariates) {
  recipe <- list(
    # Step 1: Handle missing values
    handle_missing = function(d) {
      d[complete.cases(d[c(outcome_var, treatment_var, covariates)]), ]
    },

    # Step 2: Factor treatment
    factor_treatment = function(d) {
      d[[treatment_var]] <- factor(d[[treatment_var]])
      d
    },

    # Step 3: Center covariates (for STC/MAIC)
    center_covariates = function(d, centers = NULL) {
      if (is.null(centers)) {
        centers <- sapply(d[covariates], mean, na.rm = TRUE)
      }
      for (cov in covariates) {
        d[[paste0(cov, "_centered")]] <- d[[cov]] - centers[[cov]]
      }
      attr(d, "covariate_centers") <- centers
      d
    }
  )

  class(recipe) <- c("itc_recipe", "list")
  recipe
}

# Apply recipe
prep_itc_data <- function(data, recipe) {
  result <- data
  for (step_name in names(recipe)) {
    result <- recipe[[step_name]](result)
  }
  result
}

# Unified analysis interface
run_itc_analysis <- function(
  method = c("pairwise_ma", "nma", "maic", "stc", "ml_nmr"),
  ...
) {
  method <- match.arg(method)

  # Dispatch to appropriate function
  result <- switch(method,
    pairwise_ma = run_pairwise_ma(...),
    nma = run_nma(...),
    maic = run_maic(...),
    stc = run_stc(...),
    ml_nmr = run_ml_nmr(...)
  )

  # Add common metadata
  result$method <- method
  result$timestamp <- Sys.time()
  result$session_info <- sessionInfo()

  class(result) <- c("itc_result", class(result))
  result
}

# Standard result tibble format
standardize_itc_results <- function(result) {
  tibble::tibble(
    comparison = result$comparison,
    effect_measure = result$effect_measure,
    estimate = result$estimate,
    ci_lower = result$ci_lower,
    ci_upper = result$ci_upper,
    se = result$se,
    p_value = result$p_value,
    method = result$method,
    n_studies = result$n_studies %||% NA_integer_,
    n_patients = result$n_patients %||% NA_integer_,
    heterogeneity_i2 = result$i2 %||% NA_real_,
    heterogeneity_tau2 = result$tau2 %||% NA_real_
  )
}

# Set and document seed
ANALYSIS_SEED <- 12345

# Use in all stochastic operations
set.seed(ANALYSIS_SEED)
bootstrap_result <- boot::boot(..., R = 1000)

# For parallel operations
library(doRNG)
registerDoRNG(ANALYSIS_SEED)

# Use renv for package management
renv::init()
renv::snapshot()

# Document versions in output
cat("Package versions:\n")
packageVersion("meta")
packageVersion("netmeta")
packageVersion("maicplus")

# At end of analysis
sink("session_info.txt")
sessionInfo()
sink()

# Or more detailed
writeLines(capture.output(devtools::session_info()), "session_info.txt")

validate_binary_outcome <- function(data, outcome_var) {
  vals <- data[[outcome_var]]
  if (!all(vals %in% c(0, 1, NA))) {
    stop("Binary outcome must contain only 0, 1, or NA")
  }
  if (all(vals == 0, na.rm = TRUE) || all(vals == 1, na.rm = TRUE)) {
    warning("All outcomes are identical - check data")
  }
  invisible(TRUE)
}

validate_survival_outcome <- function(data, time_var, event_var) {
  if (any(data[[time_var]] < 0, na.rm = TRUE)) {
    stop("Survival times must be non-negative")
  }
  if (!all(data[[event_var]] %in% c(0, 1, NA))) {
    stop("Event indicator must be 0, 1, or NA")
  }
  invisible(TRUE)
}

validate_agd <- function(agd, required_fields) {
  missing <- setdiff(required_fields, names(agd))
  if (length(missing) > 0) {
    stop(sprintf("Missing AgD fields: %s", paste(missing, collapse = ", ")))
  }

  # Check numeric fields are positive
  numeric_fields <- c("n_total", "n_events", "mean", "sd")
  for (field in intersect(numeric_fields, names(agd))) {
    if (any(agd[[field]] < 0, na.rm = TRUE)) {
      stop(sprintf("Field '%s' contains negative values", field))
    }
  }
  invisible(TRUE)
}

# Bad
data <- data[data$age > 65, ]

# Good
AGE_THRESHOLD <- 65
data <- data[data$age > AGE_THRESHOLD, ]

# Bad
result <- maic_anchored(weights, ipd, pseudo_ipd)

# Good
stopifnot(inherits(weights, "maicplus_estimate_weights"))
stopifnot(nrow(ipd) > 0)
result <- maic_anchored(weights, ipd, pseudo_ipd)

# Bad
bootstrap_ci <- boot::boot.ci(boot_result)

# Good
set.seed(12345)
boot_result <- boot::boot(data, statistic, R = 1000)
bootstrap_ci <- boot::boot.ci(boot_result)