## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "figures/advanced-", warning = FALSE, message = FALSE, eval = FALSE # Set to TRUE when demonstrating specific features ) ## ----load-packages, eval=TRUE------------------------------------------------- library(SCMs) library(data.table) library(ggplot2) library(parallel) library(R6) # For advanced object-oriented features set.seed(2023) ## ----performance-setup-------------------------------------------------------- # # Create larger synthetic dataset for performance testing # create_large_panel <- function(n_units = 50, n_periods = 30, n_outcomes = 3) { # # units <- c("treated", paste0("control_", 1:(n_units - 1))) # periods <- 1990:(1990 + n_periods - 1) # # panel <- expand.grid(unit = units, year = periods, stringsAsFactors = FALSE) # setDT(panel) # # # Multiple outcomes # for (i in 1:n_outcomes) { # outcome_name <- paste0("outcome_", i) # panel[, (outcome_name) := 100 + (year - 1990) * runif(1, 0.5, 2) + rnorm(.N, 0, 5)] # # # Treatment effects # panel[unit == "treated" & year >= 2005, # (outcome_name) := get(outcome_name) + rnorm(1, -10, 2)] # } # # # Multiple covariates # for (i in 1:8) { # covar_name <- paste0("covar_", i) # panel[, (covar_name) := rnorm(.N, 50, 10)] # } # # return(panel) # } # # # Create performance test dataset # large_panel <- create_large_panel(n_units = 30, n_periods = 25, n_outcomes = 2) # cat(sprintf("Created dataset with %d observations\n", nrow(large_panel))) ## ----parallel-config---------------------------------------------------------- # # Detect available cores and configure parallel processing # n_cores <- parallel::detectCores() # optimal_cores <- min(n_cores - 1, 8) # Leave one core free, cap at 8 # # cat(sprintf("System has %d cores, using %d for analysis\n", n_cores, optimal_cores)) # # # Configure parallel backend # if (require("doParallel", quietly = TRUE)) { # cl <- makeCluster(optimal_cores) # registerDoParallel(cl) # # # Ensure proper cleanup # on.exit(stopCluster(cl), add = TRUE) # } # # # Memory-aware specification processing # process_large_spec_curve <- function(dataset, max_specs_per_batch = 50) { # # # Create specification combinations # covariate_specs <- list( # basic = list( # label = "basic", # operations = list( # list(var = "covar_1", compute = "mean"), # list(var = "covar_2", compute = "mean") # ) # ), # extended = list( # label = "extended", # operations = list( # list(var = "covar_1", compute = "mean"), # list(var = "covar_2", compute = "mean"), # list(var = "covar_3", compute = "mean"), # list(var = "covar_4", compute = "mean") # ) # ), # behavioral = list( # label = "behavioral", # operations = list( # list(var = "covar_5", compute = "mean"), # list(var = "covar_6", compute = "mean") # ) # ) # ) # # # Use memory management for large specification spaces # storage_manager <- SpecResultsManager$new( # max_memory_mb = 500, # Limit memory usage # compression_level = 9 # High compression for storage # ) # # analysis_params <- list( # outcomes = c("outcome_1", "outcome_2"), # covagg = covariate_specs, # col_name_unit_name = "unit", # name_treated_unit = "treated", # treated_period = 2005, # min_period = 1990, # end_period = 2014, # col_name_period = "year", # feature_weights = c("uniform"), # outcome_models = c("none", "ridge"), # donor_sample = "all", # constants = FALSE, # constraints = list(list(name = "simplex")), # inference_type = "placebo", # inference_config = list(verbose = TRUE) # ) # # # Process analysis through the wrapper # results <- run_spec_curve_analysis( # dataset = dataset, # params = analysis_params, # cores = optimal_cores # ) # # return(results) # } ## ----custom-constraints------------------------------------------------------- # # Define custom constraint function # create_custom_constraint <- function(constraint_type = "custom", parameters = list()) { # # if (constraint_type == "elastic_net") { # # Elastic net combines L1 and L2 penalties # return(list( # name = "L1-L2", # Q1 = parameters$alpha * parameters$lambda, # L1 penalty # Q2 = (1 - parameters$alpha) * parameters$lambda # L2 penalty # )) # } # # if (constraint_type == "group_sparsity") { # # Encourage sparsity in groups of similar units # return(list( # name = "custom_group", # penalty_matrix = parameters$group_penalty_matrix, # lambda = parameters$lambda # )) # } # # if (constraint_type == "bounded_simplex") { # # Simplex with additional bounds on individual weights # return(list( # name = "bounded_simplex", # w_min = parameters$min_weight, # w_max = parameters$max_weight # )) # } # # stop(paste("Unknown constraint type:", constraint_type)) # } # # # Example custom constraints # custom_constraints <- list( # # # Elastic net with different mixing parameters # create_custom_constraint("elastic_net", # list(alpha = 0.3, lambda = 0.1)), # # create_custom_constraint("elastic_net", # list(alpha = 0.7, lambda = 0.05)), # # # Bounded simplex (prevent any single unit from dominating) # create_custom_constraint("bounded_simplex", # list(min_weight = 0.0, max_weight = 0.4)) # ) # # print("Custom constraint specifications:") # str(custom_constraints) ## ----solver-config------------------------------------------------------------ # # Intelligent solver selection based on problem characteristics # configure_optimal_solver <- function(n_donors, constraint_type, problem_size) { # # available_solvers <- CVXR::installed_solvers() # # # Commercial solvers (if available) for large problems # if (problem_size > 1000 && "GUROBI" %in% available_solvers) { # return(list( # primary = "GUROBI", # settings = list(TimeLimit = 300, MIPGap = 1e-4) # )) # } # # # Specialized solvers for different constraint types # if (constraint_type == "lasso" && "OSQP" %in% available_solvers) { # return(list( # primary = "OSQP", # settings = list(eps_abs = 1e-5, eps_rel = 1e-5) # )) # } # # if (constraint_type == "simplex" && n_donors <= 50 && "ECOS" %in% available_solvers) { # return(list( # primary = "ECOS", # settings = list(abstol = 1e-6, reltol = 1e-6) # )) # } # # # Default fallback # return(list( # primary = available_solvers[1], # settings = list() # )) # } # # # Apply solver configuration # solver_config <- configure_optimal_solver( # n_donors = 29, # constraint_type = "simplex", # problem_size = nrow(large_panel) # ) # # cat("Recommended solver configuration:\n") # str(solver_config) ## ----memory-management-------------------------------------------------------- # # Initialize advanced storage manager # storage_manager <- SpecResultsManager$new( # storage_dir = tempfile("scm_results"), # max_memory_mb = 1000, # compression_level = 6 # ) # # # Demonstrate storage and retrieval # demo_spec_results <- list( # results = data.table( # spec_id = paste0("spec_", 1:100), # treatment_effect = rnorm(100, -5, 2), # rmse = runif(100, 0.5, 2.0) # ), # abadie_inference = list( # p_values = data.table(spec_id = paste0("spec_", 1:100), p_value = runif(100)) # ) # ) # # # Store results # storage_id <- "demo_analysis_2023" # storage_manager$store_results(demo_spec_results, storage_id) # # # Retrieve specific subsets without loading all data # filtered_results <- storage_manager$get_filtered_results( # storage_id, # filter_function = function(dt) dt[treatment_effect < -3], # select_cols = c("spec_id", "treatment_effect", "rmse") # ) # # print("Filtered results sample:") # print(head(filtered_results)) # # # Get storage statistics # storage_stats <- storage_manager$get_storage_summary() # print("Storage summary:") # str(storage_stats) ## ----batch-processing--------------------------------------------------------- # # Process specifications in memory-efficient batches # process_specifications_in_batches <- function(all_specifications, batch_size = 100) { # # n_specs <- length(all_specifications) # n_batches <- ceiling(n_specs / batch_size) # # cat(sprintf("Processing %d specifications in %d batches\n", n_specs, n_batches)) # # results <- vector("list", n_batches) # # for (batch_i in 1:n_batches) { # # start_idx <- (batch_i - 1) * batch_size + 1 # end_idx <- min(batch_i * batch_size, n_specs) # # batch_specs <- all_specifications[start_idx:end_idx] # # cat(sprintf("Processing batch %d/%d (%d specifications)...\n", # batch_i, n_batches, length(batch_specs))) # # # Process batch with error handling # batch_result <- tryCatch({ # # # Your specification processing logic here # process_specification_batch(batch_specs) # # }, error = function(e) { # warning(sprintf("Batch %d failed: %s", batch_i, e$message)) # return(NULL) # }) # # results[[batch_i]] <- batch_result # # # Periodic cleanup # if (batch_i %% 5 == 0) { # gc() # Force garbage collection # } # } # # # Combine successful batches # successful_results <- results[!sapply(results, is.null)] # return(do.call(rbind, successful_results)) # } ## ----custom-plots------------------------------------------------------------- # # Create highly customized specification curve visualization # create_advanced_spec_plot <- function(spec_results, # highlight_specs = NULL, # color_by = "shap_values", # facet_by = NULL) { # # # Extract plotting data # plot_data <- spec_results$results[unit_type == "treated" & post_period == TRUE] # # # Base plot # p <- ggplot(plot_data, aes(x = reorder(full_spec_id, tau), y = tau)) # # # Color scheme based on user choice # if (color_by == "shap_values" && "shap_feature_importance" %in% names(plot_data)) { # p <- p + geom_point(aes(color = shap_feature_importance), alpha = 0.7, size = 1.5) # p <- p + scale_color_gradient2( # low = "blue", mid = "white", high = "red", # name = "SHAP\nImportance" # ) # } else if (color_by == "rmse") { # p <- p + geom_point(aes(color = rmse), alpha = 0.7, size = 1.5) # p <- p + scale_color_gradient(low = "darkgreen", high = "darkred", name = "RMSE") # } else { # p <- p + geom_point(alpha = 0.6, size = 1.2, color = "steelblue") # } # # # Highlight specific specifications # if (!is.null(highlight_specs)) { # highlight_data <- plot_data[full_spec_id %in% highlight_specs] # p <- p + geom_point(data = highlight_data, color = "red", size = 3, shape = 1, stroke = 2) # } # # # Add reference lines # p <- p + geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) # p <- p + geom_hline(yintercept = median(plot_data$tau, na.rm = TRUE), # color = "darkred", linetype = "solid", alpha = 0.7) # # # Faceting if requested # if (!is.null(facet_by) && facet_by %in% names(plot_data)) { # p <- p + facet_wrap(as.formula(paste("~", facet_by)), scales = "free_x") # } # # # Styling # p <- p + theme_minimal() + # theme( # axis.text.x = element_blank(), # axis.ticks.x = element_blank(), # panel.grid.minor = element_blank(), # plot.title = element_text(size = 14, hjust = 0.5), # plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray50") # ) + # labs( # x = "Specifications (ordered by effect size)", # y = "Treatment Effect Estimate", # title = "Advanced Specification Curve Analysis", # subtitle = paste("Distribution of treatment effects across", # nrow(plot_data), "specifications") # ) # # return(p) # } # # # Example usage (with mock data) # mock_results <- list( # results = data.table( # full_spec_id = paste0("spec_", 1:200), # unit_type = "treated", # post_period = TRUE, # tau = rnorm(200, -8, 3), # rmse = runif(200, 0.8, 2.5), # outcome_model = sample(c("none", "ridge", "lasso"), 200, replace = TRUE), # const = sample(c("simplex", "lasso", "ridge"), 200, replace = TRUE) # ) # ) # # # Create custom plot # advanced_plot <- create_advanced_spec_plot( # mock_results, # color_by = "rmse", # facet_by = "outcome_model" # ) # # print(advanced_plot) ## ----shap-analysis------------------------------------------------------------ # # Advanced SHAP analysis functions # analyze_shap_patterns <- function(shap_data, specification_data) { # # if (is.null(shap_data) || nrow(shap_data) == 0) { # warning("No SHAP data available for analysis") # return(NULL) # } # # # Merge SHAP values with specification results # merged_data <- merge(shap_data, specification_data, by = "full_spec_id") # # # Calculate feature importance statistics # feature_importance <- merged_data[, .( # mean_abs_shap = mean(abs(shap_value), na.rm = TRUE), # sd_shap = sd(shap_value, na.rm = TRUE), # positive_proportion = mean(shap_value > 0, na.rm = TRUE), # correlation_with_effect = cor(shap_value, tau, use = "complete.obs") # ), by = feature_name][order(-mean_abs_shap)] # # return(list( # feature_importance = feature_importance, # merged_data = merged_data # )) # } # # # Create SHAP interaction plots # plot_shap_interactions <- function(shap_analysis_result) { # # if (is.null(shap_analysis_result)) return(NULL) # # merged_data <- shap_analysis_result$merged_data # # # Create interaction heatmap # interaction_plot <- ggplot(merged_data, # aes(x = feature_name, y = reorder(full_spec_id, tau))) + # geom_tile(aes(fill = shap_value)) + # scale_fill_gradient2(low = "blue", mid = "white", high = "red", # name = "SHAP\nValue") + # theme_minimal() + # theme(axis.text.y = element_blank(), # axis.ticks.y = element_blank(), # axis.text.x = element_text(angle = 45, hjust = 1)) + # labs(x = "Features", y = "Specifications (ordered by effect size)", # title = "SHAP Value Heatmap Across Specifications") # # return(interaction_plot) # } ## ----custom-inference--------------------------------------------------------- # # Implement custom inference procedure # custom_bootstrap_inference <- function(scdata_obj, n_bootstrap = 1000, # confidence_level = 0.95) { # # # Extract data components # Y_treated_pre <- scdata_obj$Y.pre # Y_donors_pre <- scdata_obj$Y.donors # # n_pre_periods <- length(Y_treated_pre) # n_donors <- ncol(Y_donors_pre) # # # Bootstrap procedure # bootstrap_effects <- replicate(n_bootstrap, { # # # Resample time periods with replacement # bootstrap_periods <- sample(1:n_pre_periods, n_pre_periods, replace = TRUE) # # # Create bootstrap data # Y_boot_treated <- Y_treated_pre[bootstrap_periods] # Y_boot_donors <- Y_donors_pre[bootstrap_periods, , drop = FALSE] # # # Estimate synthetic control on bootstrap sample # bootstrap_scdata <- scdata_obj # bootstrap_scdata$Y.pre <- Y_boot_treated # bootstrap_scdata$Y.donors <- Y_boot_donors # # # Estimate effect (simplified for demonstration) # tryCatch({ # boot_result <- scest(bootstrap_scdata, w.constr = list(name = "simplex")) # boot_effect <- mean(boot_result$est.results$Y.post.fit - scdata_obj$Y.post) # return(boot_effect) # }, error = function(e) { # return(NA) # }) # }) # # # Calculate confidence intervals # bootstrap_effects <- bootstrap_effects[!is.na(bootstrap_effects)] # # alpha <- 1 - confidence_level # ci_lower <- quantile(bootstrap_effects, alpha/2, na.rm = TRUE) # ci_upper <- quantile(bootstrap_effects, 1 - alpha/2, na.rm = TRUE) # # return(list( # bootstrap_distribution = bootstrap_effects, # confidence_interval = c(ci_lower, ci_upper), # p_value = 2 * min(mean(bootstrap_effects >= 0, na.rm = TRUE), # mean(bootstrap_effects <= 0, na.rm = TRUE)) # )) # } ## ----feature-engineering------------------------------------------------------ # # Advanced feature engineering for synthetic control # engineer_advanced_features <- function(data, outcome_var, time_var, id_var, # feature_specs = list()) { # # setDT(data) # # # Create polynomial time trends # if ("polynomial_trends" %in% names(feature_specs)) { # degree <- feature_specs$polynomial_trends$degree # data[, paste0("time_poly_", 1:degree) := lapply(1:degree, function(d) { # (get(time_var) - min(get(time_var)))^d # })] # } # # # Create interaction terms between covariates # if ("interactions" %in% names(feature_specs)) { # interaction_pairs <- feature_specs$interactions$pairs # for (pair in interaction_pairs) { # new_var <- paste(pair, collapse = "_x_") # data[, (new_var) := get(pair[1]) * get(pair[2])] # } # } # # # Create rolling averages of outcomes # if ("rolling_outcomes" %in% names(feature_specs)) { # window_sizes <- feature_specs$rolling_outcomes$windows # for (window in window_sizes) { # new_var <- paste0(outcome_var, "_roll_", window) # data[, (new_var) := frollmean(get(outcome_var), window), by = get(id_var)] # } # } # # # Create volatility measures # if ("volatility" %in% names(feature_specs)) { # window <- feature_specs$volatility$window # vol_var <- paste0(outcome_var, "_volatility") # data[, (vol_var) := frollapply(get(outcome_var), window, sd, na.rm = TRUE), # by = get(id_var)] # } # # return(data) # } # # # Example feature engineering specification # advanced_features <- list( # polynomial_trends = list(degree = 3), # interactions = list(pairs = list(c("income", "price"), c("age", "education"))), # rolling_outcomes = list(windows = c(3, 5)), # volatility = list(window = 5) # ) # # # Apply feature engineering (example with mock data) # # engineered_data <- engineer_advanced_features( # # large_panel, "outcome_1", "year", "unit", advanced_features # # ) ## ----xgboost-integration, eval=FALSE------------------------------------------ # # Advanced SHAP analysis with XGBoost # run_advanced_shap_analysis <- function(spec_results, xgboost_config = NULL) { # # if (is.null(xgboost_config)) { # xgboost_config <- create_xgboost_config( # dataset_name = "analysis", # treated_unit_name = "treated" # ) # } # # # Run XGBoost SHAP analysis on long-format specification results # shap_results <- run_xgboost_shap_analysis( # long_data = spec_results$results, # config = xgboost_config # ) # # return(shap_results) # } ## ----database-integration, eval=FALSE----------------------------------------- # # Integration with database backends for very large datasets # setup_database_backend <- function(connection_string, table_name) { # # require(DBI) # require(RSQLite) # or appropriate database driver # # # Establish connection # con <- dbConnect(RSQLite::SQLite(), connection_string) # # # Create optimized table structure # dbExecute(con, sprintf(" # CREATE TABLE IF NOT EXISTS %s ( # unit_id TEXT, # time_period INTEGER, # outcome_value REAL, # covariate_1 REAL, # covariate_2 REAL, # /* ... additional columns ... */ # INDEX(unit_id, time_period) # ) # ", table_name)) # # return(con) # } # # # Query data efficiently for specification curve analysis # query_specification_data <- function(db_connection, specification_params) { # # # Build optimized query based on specification # query <- sprintf(" # SELECT unit_id, time_period, outcome_value, %s # FROM panel_data # WHERE time_period BETWEEN %d AND %d # AND unit_id IN (%s) # ", # paste(specification_params$covariates, collapse = ", "), # specification_params$min_period, # specification_params$max_period, # paste0("'", specification_params$units, "'", collapse = ", ") # ) # # # Execute query and return data.table # result <- dbGetQuery(db_connection, query) # return(as.data.table(result)) # } ## ----memory-profiling--------------------------------------------------------- # # Monitor memory usage during large analyses # monitor_memory_usage <- function(analysis_function, ...) { # # # Baseline memory # baseline_memory <- pryr::mem_used() # # cat(sprintf("Baseline memory usage: %s\n", # format(baseline_memory, units = "auto"))) # # # Run analysis with memory monitoring # start_time <- Sys.time() # # result <- tryCatch({ # analysis_function(...) # }, error = function(e) { # cat(sprintf("Analysis failed: %s\n", e$message)) # return(NULL) # }) # # end_time <- Sys.time() # peak_memory <- pryr::mem_used() # # # Report results # cat(sprintf("Analysis completed in: %s\n", # format(end_time - start_time))) # cat(sprintf("Peak memory usage: %s\n", # format(peak_memory, units = "auto"))) # cat(sprintf("Additional memory used: %s\n", # format(peak_memory - baseline_memory, units = "auto"))) # # # Cleanup # gc() # # return(result) # } ## ----optimization-tips-------------------------------------------------------- # # Provide optimization recommendations based on data characteristics # recommend_optimizations <- function(n_units, n_periods, n_specifications, # available_memory_gb = 8) { # # # Estimate memory requirements # estimated_memory_gb <- (n_units * n_periods * n_specifications * 8) / (1024^3) # # recommendations <- list() # # # Memory recommendations # if (estimated_memory_gb > available_memory_gb) { # recommendations$memory <- c( # "Consider batch processing", # "Use disk-based storage manager", # "Increase compression level", # "Reduce specification space if possible" # ) # } # # # Computational recommendations # if (n_specifications > 500) { # recommendations$computation <- c( # "Enable parallel processing", # "Consider distributed computing", # "Use faster solvers if available" # ) # } # # # Data structure recommendations # if (n_units > 100) { # recommendations$data_structure <- c( # "Use data.table for all operations", # "Consider database backend", # "Implement sparse matrix representations where possible" # ) # } # # return(recommendations) # } # # # Example usage # opts <- recommend_optimizations( # n_units = 200, # n_periods = 40, # n_specifications = 1000, # available_memory_gb = 16 # ) # # cat("Optimization recommendations:\n") # str(opts)