## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)
library(WMAP)

## ----install, eval=FALSE------------------------------------------------------
# # Install from CRAN
# install.packages("WMAP")
# 
# # Load package
# library(WMAP)

## ----load-data----------------------------------------------------------------
# Load demonstration dataset
data(demo)

# The demo data contains:
# S - Study/site indicators (7 hospitals)
# Z - Group indicators (2 groups: IDC vs ILC)
# X - Patient covariates (clinical variables)
# Y - Outcomes (8 gene expression measurements)
# naturalGroupProp - Known population group proportions

## ----basic-analysis-----------------------------------------------------------
# Estimate causal effects
result <- causal.estimate(
  S = S,                           # Study indicators
  Z = Z,                           # Group indicators
  X = X,                           # Covariates
  Y = Y,                           # Outcomes
  B = 20,                          # Bootstrap samples (use 100+ in practice)
  method = "IC",                   # Weighting method
  naturalGroupProp = naturalGroupProp,
  seed = 123                       # For reproducibility
)

## ----results------------------------------------------------------------------
# Print summary
summary(result)

# Visualize bootstrap ESS distribution
plot(result)

## ----methods, eval=FALSE------------------------------------------------------
# # Compare methods
# result_ic <- causal.estimate(S, Z, X, Y, B = 100, method = "IC", naturalGroupProp)
# result_igo <- causal.estimate(S, Z, X, Y, B = 100, method = "IGO", naturalGroupProp)
# result_flexor <- causal.estimate(S, Z, X, Y, B = 100, method = "FLEXOR", naturalGroupProp)
# 
# # Compare ESS
# result_ic$percentESS
# result_igo$percentESS
# result_flexor$percentESS

## ----outcome-model, eval=FALSE------------------------------------------------
# # Default: use observed outcomes directly
# result <- causal.estimate(S, Z, X, Y, B = 100, method = "IC",
#                           naturalGroupProp = naturalGroupProp)
# 
# # RF predictions substituted for observed Y
# result <- causal.estimate(S, Z, X, Y, B = 100, method = "IC",
#                           naturalGroupProp = naturalGroupProp,
#                           outcome_model = TRUE)

## ----crossfit, eval=FALSE-----------------------------------------------------
# # Enable 5-fold cross-fitting
# result <- causal.estimate(S, Z, X, Y, B = 100, method = "IC",
#                           naturalGroupProp = naturalGroupProp,
#                           outcome_model = TRUE,
#                           crossfit_folds = 5)

## ----bootstrap-size, eval=FALSE-----------------------------------------------
# # Minimum for stable CIs
# result <- causal.estimate(S, Z, X, Y, B = 50, method = "IC", naturalGroupProp)
# 
# # Recommended for publication
# result <- causal.estimate(S, Z, X, Y, B = 200, method = "IC", naturalGroupProp)
# 
# # High precision
# result <- causal.estimate(S, Z, X, Y, B = 500, method = "IC", naturalGroupProp)

## ----no-bootstrap, eval=FALSE-------------------------------------------------
# # Point estimates only (no confidence intervals)
# result <- causal.estimate(S, Z, X, Y, B = 0, method = "IC", naturalGroupProp)
# summary(result)  # Shows point estimates without CIs

## ----parallel, eval=FALSE-----------------------------------------------------
# # Install required packages first
# install.packages(c("future", "future.apply"))
# 
# # Define n_cores based on your system (e.g., 4)
# result <- causal.estimate(
#   S, Z, X, Y,
#   B = 100,
#   method = "IC",
#   naturalGroupProp = naturalGroupProp,
#   parallel = TRUE,
#   n_cores = 4
# )

## ----low-ess-demo, eval=FALSE-------------------------------------------------
# # Example warning:
# # Warning message:
# # Low effective sample size: 3.1% (threshold: 5%).

## ----extreme-weights-demo, eval=FALSE-----------------------------------------
# # Example warning:
# # Warning message:
# # Extreme weights: max/mean ratio = 45.2 (threshold: 20).

## ----positivity-demo, eval=FALSE----------------------------------------------
# # Example warning:
# # Warning message:
# # Sparse study-group cells: 2 cell(s) with < 5 observations.
# 
# # Example warning:
# # Warning message:
# # Near-zero propensity scores: 6.3% of subjects truncated.

## ----standard-workflow, eval=FALSE--------------------------------------------
# # 1. Load and check data
# data(demo)
# table(S, Z)  # Check all cells have observations
# sum(is.na(cbind(S, Z, X, Y)))  # Check for missing values
# 
# # 2. Run analysis
# result <- causal.estimate(
#   S, Z, X, Y,
#   B = 100,
#   method = "FLEXOR",
#   naturalGroupProp = naturalGroupProp,
#   seed = 123
# )
# 
# # 3. Check diagnostics
# summary(result)
# plot(result)

## ----comparison-workflow, eval=FALSE------------------------------------------
# # Compare all three methods
# set.seed(123)
# 
# results_ic <- causal.estimate(S, Z, X, Y, B = 100, method = "IC", naturalGroupProp)
# results_igo <- causal.estimate(S, Z, X, Y, B = 100, method = "IGO", naturalGroupProp)
# results_flexor <- causal.estimate(S, Z, X, Y, B = 100, method = "FLEXOR", naturalGroupProp)
# 
# # Compare ESS
# cat("IC ESS:", results_ic$percentESS, "%\n")
# cat("IGO ESS:", results_igo$percentESS, "%\n")
# cat("FLEXOR ESS:", results_flexor$percentESS, "%\n")
# 
# # Compare estimates
# summary(results_ic)
# summary(results_igo)
# summary(results_flexor)

## ----sensitivity-workflow, eval=FALSE-----------------------------------------
# # Test sensitivity to bootstrap sample size
# result_b50 <- causal.estimate(S, Z, X, Y, B = 50, method = "IGO", naturalGroupProp, seed = 123)
# result_b100 <- causal.estimate(S, Z, X, Y, B = 100, method = "IGO", naturalGroupProp, seed = 123)
# result_b200 <- causal.estimate(S, Z, X, Y, B = 200, method = "IGO", naturalGroupProp, seed = 123)
# 
# # Compare CI widths
# summary(result_b50)
# summary(result_b100)
# summary(result_b200)

## ----weights-workflow, eval=FALSE---------------------------------------------
# # Get weights
# weights <- balancing.weights(
#   S, Z, X,
#   method = "IC",
#   naturalGroupProp = naturalGroupProp
# )
# 
# # Examine weight distribution
# summary(weights)
# 
# # Access weight vector for custom analysis
# head(weights$wt.v)

## ----reproducibility, eval=FALSE----------------------------------------------
# result1 <- causal.estimate(S, Z, X, Y, B = 100, method = "IC", naturalGroupProp, seed = 123)
# result2 <- causal.estimate(S, Z, X, Y, B = 100, method = "IC", naturalGroupProp, seed = 123)

## ----data-quality, eval=FALSE-------------------------------------------------
# # Check for missing values
# sum(is.na(cbind(S, Z, X, Y)))
# 
# # Ensure all study-group combinations have observations
# table(S, Z)

## ----help, eval=FALSE---------------------------------------------------------
# # Function documentation
# ?causal.estimate
# ?balancing.weights
# 
# # Package citation
# citation("WMAP")