## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "figures/workflow-",
  warning = FALSE,
  message = FALSE,
  eval = FALSE  # Set to TRUE when running with actual data
)

## ----load-packages, eval=TRUE-------------------------------------------------
library(SCMs)
library(data.table)
library(ggplot2)
library(dplyr)

set.seed(42)  # For reproducible results

## ----load-data----------------------------------------------------------------
# # In practice, you would load your data here
# # data(tobacco_data)  # Example dataset
# # For this vignette, we'll create realistic synthetic data
# 
# # Create realistic tobacco consumption data
# create_tobacco_data <- function() {
# 
#   states <- c("California", "Texas", "New York", "Florida", "Pennsylvania",
#               "Illinois", "Ohio", "Georgia", "North Carolina", "Michigan",
#               "New Jersey", "Virginia", "Washington", "Arizona", "Massachusetts")
# 
#   years <- 1980:2000
# 
#   # Create panel structure
#   panel <- expand.grid(state = states, year = years, stringsAsFactors = FALSE)
#   setDT(panel)
# 
#   # Add realistic covariates
#   panel[, `:=`(
#     # Per capita cigarette consumption (packs per year)
#     cigsale = 120 - (year - 1980) * 1.5 + rnorm(.N, 0, 8),
# 
#     # Retail price per pack (dollars)
#     retprice = 0.60 + (year - 1980) * 0.05 + rnorm(.N, 0, 0.05),
# 
#     # Per capita income (thousands of dollars)
#     income = 15 + (year - 1980) * 0.8 + rnorm(.N, 0, 1.5),
# 
#     # Population 15 years and older (millions)
#     age15to24 = 2.5 + rnorm(.N, 0, 0.3),
# 
#     # Beer consumption per capita
#     beer = 22 + rnorm(.N, 0, 2)
#   )]
# 
#   # Add California treatment effect after 1988
#   panel[state == "California" & year >= 1989,
#         cigsale := cigsale - (year - 1988) * 8]
# 
#   # Ensure realistic bounds
#   panel[cigsale < 0, cigsale := 0]
#   panel[retprice < 0.30, retprice := 0.30]
#   panel[income < 8, income := 8]
# 
#   return(panel)
# }
# 
# # Generate the data
# tobacco_data <- create_tobacco_data()
# 
# # Examine the data structure
# str(tobacco_data)
# head(tobacco_data)
# 
# # Summary statistics
# summary(tobacco_data)

## ----exploratory-analysis-----------------------------------------------------
# # Plot cigarette consumption trends
# ggplot(tobacco_data, aes(x = year, y = cigsale, color = state)) +
#   geom_line(alpha = 0.7) +
#   geom_line(data = tobacco_data[state == "California"],
#             color = "red", size = 1.2) +
#   geom_vline(xintercept = 1988, linetype = "dashed", color = "red", alpha = 0.7) +
#   labs(
#     title = "Per Capita Cigarette Consumption by State",
#     subtitle = "California highlighted in red, treatment begins 1988",
#     x = "Year",
#     y = "Cigarette Packs Per Capita",
#     color = "State"
#   ) +
#   theme_minimal() +
#   theme(legend.position = "none")
# 
# # Pre-treatment comparison
# pre_treatment_avg <- tobacco_data[year <= 1988,
#   .(avg_consumption = mean(cigsale, na.rm = TRUE)),
#   by = state
# ][order(-avg_consumption)]
# 
# print("Pre-treatment average cigarette consumption (1980-1988):")
# print(pre_treatment_avg)
# 
# # Check for missing data
# tobacco_data[, lapply(.SD, function(x) sum(is.na(x)))]

## ----covariate-specs----------------------------------------------------------
# # Define multiple covariate specifications to test
# covariate_specifications <- list(
# 
#   # Specification 1: Economic factors
#   economic = list(
#     label = "Economic Variables",
#     operations = list(
#       list(var = "retprice", compute = "mean"),
#       list(var = "income", compute = "mean")
#     )
#   ),
# 
#   # Specification 2: Full demographic and economic
#   full_demographics = list(
#     label = "Full Demographics + Economics",
#     operations = list(
#       list(var = "retprice", compute = "mean"),
#       list(var = "income", compute = "mean"),
#       list(var = "age15to24", compute = "mean"),
#       list(var = "beer", compute = "mean")
#     )
#   ),
# 
#   # Specification 3: Outcome lags
#   with_lags = list(
#     label = "With Outcome Lags",
#     operations = list(
#       list(var = "outcome_var", partition_periods = list(type = "by_period")),
#       list(var = "retprice", compute = "mean"),
#       list(var = "income", compute = "mean")
#     )
#   ),
# 
#   # Specification 4: Price focus
#   price_focused = list(
#     label = "Price-Focused",
#     operations = list(
#       list(var = "retprice", compute = "mean"),
#       list(var = "retprice", partition_periods = list(type = "by_period"))
#     )
#   ),
# 
#   # Specification 5: Behavioral proxies
#   behavioral = list(
#     label = "Behavioral Proxies",
#     operations = list(
#       list(var = "beer", compute = "mean"),
#       list(var = "age15to24", compute = "mean")
#     )
#   ),
# 
#   # Specification 6: Minimal (for robustness)
#   minimal = list(
#     label = "Outcome Only",
#     operations = list(
#       list(var = "outcome_var", partition_periods = list(type = "by_period"))
#     )
#   )
# )
# 
# # Display specification details
# cat("Covariate Specifications:\n")
# for (i in seq_along(covariate_specifications)) {
#   spec <- covariate_specifications[[i]]
#   cat(sprintf("\n%d. %s:\n", i, spec$label))
#   cat(sprintf("   Operations: %d\n", length(spec$operations)))
# }

## ----config-options-----------------------------------------------------------
# # Define the full specification space to explore
# analysis_config <- list(
# 
#   # Feature weighting methods
#   feature_weights = c("uniform", "optimize"),
# 
#   # Outcome modeling approaches
#   outcome_models = c("none", "augsynth", "ridge", "lasso"),
# 
#   # Weight constraints
#   constraints = list(
#     list(name = "simplex"),                    # Standard SC
#     list(name = "lasso", Q = 0.01),           # L1 regularized
#     list(name = "lasso", Q = 0.1),            # More L1 regularization
#     list(name = "ridge", Q = 0.01),           # L2 regularized
#     list(name = "ols")                        # Unconstrained
#   ),
# 
#   # Donor selection methods
#   donor_sample = c("all")  # Use all available control states
# )
# 
# # Print configuration summary
# cat("Analysis Configuration:\n")
# cat(sprintf("- Feature weighting methods: %d\n", length(analysis_config$feature_weights)))
# cat(sprintf("- Outcome models: %d\n", length(analysis_config$outcome_models)))
# cat(sprintf("- Weight constraints: %d\n", length(analysis_config$constraints)))
# cat(sprintf("- Covariate specifications: %d\n", length(covariate_specifications)))
# 
# # Calculate total specifications
# total_specs <- length(covariate_specifications) *
#                length(analysis_config$feature_weights) *
#                length(analysis_config$outcome_models) *
#                length(analysis_config$constraints)
# 
# cat(sprintf("\nTotal specifications to estimate: %d\n", total_specs))

## ----run-analysis-------------------------------------------------------------
# # Run the comprehensive specification curve analysis
# cat("Running specification curve analysis...\n")
# cat("This may take several minutes depending on the number of specifications.\n\n")
# 
# analysis_params <- list(
#   outcomes = "cigsale",
#   covagg = covariate_specifications,
#   col_name_unit_name = "state",
#   name_treated_unit = "California",
#   col_name_period = "year",
#   treated_period = 1989,  # Treatment began in 1989
#   min_period = 1980,
#   end_period = 2000,
#   feature_weights = analysis_config$feature_weights,
#   outcome_models = analysis_config$outcome_models,
#   constraints = analysis_config$constraints,
#   donor_sample = analysis_config$donor_sample,
#   constants = FALSE,
#   inference_type = "all",  # Both placebo and bootstrap
#   inference_config = list(
#     bootstrap_n_replications = 500,
#     bootstrap_cores = 1
#   )
# )
# 
# tobacco_results <- run_spec_curve_analysis(
#   dataset = tobacco_data,
#   params = analysis_params,
#   cores = 2  # Adjust based on your system
# )
# 
# # Check if analysis completed successfully
# if (!is.null(tobacco_results)) {
#   cat("Analysis completed successfully!\n")
#   cat(sprintf("Number of successful specifications: %d\n",
#               nrow(tobacco_results$results)))
# } else {
#   stop("Analysis failed to complete")
# }

## ----results-summary----------------------------------------------------------
# # Print overall summary
# print(tobacco_results)
# summary(tobacco_results)
# 
# # Extract key statistics
# results_dt <- tobacco_results$results
# treated_effects <- results_dt[unit_type == "treated" & post_period == TRUE]
# 
# cat("Treatment Effect Summary:\n")
# cat(sprintf("  Mean effect: %.2f packs per capita\n", mean(treated_effects$tau, na.rm = TRUE)))
# cat(sprintf("  Median effect: %.2f packs per capita\n", median(treated_effects$tau, na.rm = TRUE)))
# cat(sprintf("  Standard deviation: %.2f\n", sd(treated_effects$tau, na.rm = TRUE)))
# cat(sprintf("  Range: [%.2f, %.2f]\n",
#             min(treated_effects$tau, na.rm = TRUE),
#             max(treated_effects$tau, na.rm = TRUE)))
# 
# # Calculate percentage reduction
# baseline_consumption <- mean(tobacco_data[state == "California" & year <= 1988, cigsale])
# mean_effect <- mean(treated_effects$tau, na.rm = TRUE)
# percent_reduction <- abs(mean_effect) / baseline_consumption * 100
# 
# cat(sprintf("  Estimated effect size: %.1f%% reduction in consumption\n", percent_reduction))

## ----plot-spec-curve----------------------------------------------------------
# # Main specification curve plot
# spec_plot <- plot(tobacco_results)
# print(spec_plot)
# 
# # Customized plot with additional annotations
# custom_plot <- plot_spec_curve(
#   tobacco_results,
#   title = "California Tobacco Control Program: Specification Curve Analysis",
#   subtitle = sprintf("Effect on per-capita cigarette consumption (%d specifications)",
#                      nrow(treated_effects)),
#   x_axis_label = "Specification (ranked by effect size)",
#   y_axis_label = "Treatment Effect (packs per capita per year)",
#   confidence_level = 0.95
# )
# 
# print(custom_plot)

## ----inference-results--------------------------------------------------------
# # Examine inference results
# if ("abadie_inference" %in% names(tobacco_results)) {
# 
#   abadie_results <- tobacco_results$abadie_inference
# 
#   if ("p_values_rmse_ratio" %in% names(abadie_results)) {
#     p_values <- abadie_results$p_values_rmse_ratio
# 
#     cat("Abadie Placebo Test Results:\n")
#     cat(sprintf("  Specifications with p < 0.05: %d (%.1f%%)\n",
#                 sum(p_values$p_value < 0.05, na.rm = TRUE),
#                 mean(p_values$p_value < 0.05, na.rm = TRUE) * 100))
#     cat(sprintf("  Specifications with p < 0.10: %d (%.1f%%)\n",
#                 sum(p_values$p_value < 0.10, na.rm = TRUE),
#                 mean(p_values$p_value < 0.10, na.rm = TRUE) * 100))
#     cat(sprintf("  Median p-value: %.3f\n", median(p_values$p_value, na.rm = TRUE)))
#   }
# }
# 
# if ("bootstrap_inference" %in% names(tobacco_results)) {
# 
#   bootstrap_results <- tobacco_results$bootstrap_inference
# 
#   if ("p_values" %in% names(bootstrap_results)) {
#     boot_p_values <- bootstrap_results$p_values
# 
#     cat("\nBootstrap Test Results:\n")
#     cat(sprintf("  Specifications with p < 0.05: %d (%.1f%%)\n",
#                 sum(boot_p_values$p_value_two_tailed < 0.05, na.rm = TRUE),
#                 mean(boot_p_values$p_value_two_tailed < 0.05, na.rm = TRUE) * 100))
#     cat(sprintf("  Median p-value: %.3f\n",
#                 median(boot_p_values$p_value_two_tailed, na.rm = TRUE)))
#   }
# }

## ----effect-analysis----------------------------------------------------------
# # Analyze effect patterns by specification characteristics
# results_with_specs <- treated_effects
# 
# # Add specification characteristics for analysis
# # (In practice, this would be done automatically by the package)
# results_with_specs[, `:=`(
#   is_augmented = grepl("augsynth|ridge|lasso", outcome_model),
#   is_regularized = grepl("lasso|ridge", const),
#   uses_optimization = fw == "optimize"
# )]
# 
# # Effect by outcome model
# outcome_model_effects <- results_with_specs[,
#   .(mean_effect = mean(tau, na.rm = TRUE),
#     sd_effect = sd(tau, na.rm = TRUE),
#     n_specs = .N),
#   by = outcome_model
# ][order(mean_effect)]
# 
# cat("Treatment effects by outcome model:\n")
# print(outcome_model_effects)
# 
# # Effect by constraint type
# constraint_effects <- results_with_specs[,
#   .(mean_effect = mean(tau, na.rm = TRUE),
#     sd_effect = sd(tau, na.rm = TRUE),
#     n_specs = .N),
#   by = const
# ][order(mean_effect)]
# 
# cat("\nTreatment effects by constraint type:\n")
# print(constraint_effects)

## ----robustness---------------------------------------------------------------
# # Assess robustness using multiple criteria
# 
# # 1. Sign consistency
# negative_effects <- mean(treated_effects$tau < 0, na.rm = TRUE)
# cat(sprintf("Proportion of specifications with negative effects: %.1f%%\n",
#             negative_effects * 100))
# 
# # 2. Statistical significance across specifications
# if (exists("p_values")) {
#   significant_specs <- mean(p_values$p_value < 0.05, na.rm = TRUE)
#   cat(sprintf("Proportion of statistically significant specifications: %.1f%%\n",
#               significant_specs * 100))
# }
# 
# # 3. Effect size stability
# effect_cv <- sd(treated_effects$tau, na.rm = TRUE) / abs(mean(treated_effects$tau, na.rm = TRUE))
# cat(sprintf("Coefficient of variation in effect sizes: %.2f\n", effect_cv))
# 
# # 4. Outlier analysis
# q1 <- quantile(treated_effects$tau, 0.25, na.rm = TRUE)
# q3 <- quantile(treated_effects$tau, 0.75, na.rm = TRUE)
# iqr <- q3 - q1
# outlier_threshold <- 1.5 * iqr
# 
# outliers_low <- sum(treated_effects$tau < (q1 - outlier_threshold), na.rm = TRUE)
# outliers_high <- sum(treated_effects$tau > (q3 + outlier_threshold), na.rm = TRUE)
# outlier_proportion <- (outliers_low + outliers_high) / nrow(treated_effects)
# 
# cat(sprintf("Proportion of outlier specifications: %.1f%%\n", outlier_proportion * 100))
# 
# # 5. Robustness conclusion
# if (negative_effects > 0.8 && effect_cv < 0.5) {
#   cat("\n✓ CONCLUSION: Results appear robust across specifications\n")
# } else if (negative_effects > 0.6 && effect_cv < 0.8) {
#   cat("\n⚠ CONCLUSION: Results show moderate robustness\n")
# } else {
#   cat("\n⚠ CONCLUSION: Results show limited robustness - interpret with caution\n")
# }

## ----shap-analysis, eval=FALSE------------------------------------------------
# # Note: This section would require SHAP values to be computed
# # This is typically done as part of the specification curve analysis
# 
# # if ("shap_values" %in% names(tobacco_results)) {
# #
# #   # Analyze which features drive variation in treatment effects
# #   shap_summary <- analyze_shap_importance(tobacco_results$shap_values)
# #
# #   # Plot SHAP value distributions
# #   plot_shap_distributions(tobacco_results$shap_values)
# #
# #   # Feature interaction analysis
# #   plot_shap_interactions(tobacco_results$shap_values)
# # }