Bootstrap and Power Analysis

library(tna)

# Generate some data
sequences <- simulate_sequences(
  n_sequences = 150,
  seq_length = 25,
  n_states = 5,
  seed = 42
)

# Fit original model
original_model <- fit_network_model(sequences, "tna")

# Run bootstrap
bootstrap_results <- run_bootstrap_simulation(
  original_data = sequences,
  n_bootstrap = 100,
  model_type = "tna",
  sample_size = NULL,  # Use original size
  parallel = TRUE,
  seed = 42
)

# Evaluate bootstrap performance
evaluation <- run_bootstrap_iteration(
  bootstrap_results = bootstrap_results,
  original_model = original_model
)

# View summary
print(evaluation)

# Access specific metrics
evaluation$edge_stability
evaluation$centrality_stability
evaluation$confidence_intervals

# Define stable transitions to test
stable_transitions <- list(
  c("Plan", "Monitor"),
  c("Monitor", "Evaluate")
)

# Bootstrap with stability testing
bootstrap_results <- run_bootstrap_simulation(
  original_data = sequences,
  n_bootstrap = 100,
  model_type = "tna",
  stable_transitions = stable_transitions,
  seed = 42
)

# Create parameter grid
power_grid <- generate_param_grid(
  param_ranges = list(
    n_sequences = c(25, 300),
    seq_length = c(15, 40),
    n_states = c(4, 8)
  ),
  n = 50,
  method = "lhs"
)

nrow(power_grid)  # Number of conditions
head(power_grid)

# Run power analysis
power_results <- run_grid_simulation(
  param_grid = power_grid,
  n_runs_per_setting = 50,  # Replications per condition
  model_type = "tna",
  parallel = TRUE,
  seed = 42
)

# Comprehensive analysis
analysis <- summarize_grid_results(power_results)

# Results by sample size
analysis$by_n_sequences

# Results by sequence length
analysis$by_seq_length

# Results by number of states
analysis$by_n_states

# Interaction effects
analysis$interactions

# Find minimum n_sequences for target correlation
target_correlation <- 0.90

summary_by_n <- summarize_simulation(
  results = power_results,
  by = "n_sequences",
  metrics = c("mean", "ci_lower")
)

# Find first sample size meeting target
min_n <- summary_by_n$n_sequences[
  which(summary_by_n$correlation_mean >= target_correlation)[1]
]

cat("Minimum sample size for r >=", target_correlation, ":", min_n)

library(ggplot2)

ggplot(power_df, aes(x = n_sequences, y = mean_r)) +
geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
  geom_hline(yintercept = 0.90, linetype = "dashed", color = "red") +
  labs(
    x = "Number of Sequences",
    y = "Mean Correlation with True Network",
    title = "Power Analysis: Network Recovery by Sample Size"
  ) +
  theme_minimal()

# For each sample size, calculate edge recovery
recovery_results <- lapply(sample_sizes, function(n) {
  recoveries <- lapply(1:50, function(rep) {
    seqs <- simulate_sequences(
      trans_matrix = true_probs$transition_matrix,
      init_probs = true_probs$initial_probs,
      n_sequences = n,
      seq_length = 25,
      seed = rep * 1000 + n
    )

    model <- fit_network_model(seqs, "tna")

    # Create pseudo-model from true matrix
    true_model <- list(weights = true_probs$transition_matrix)

    recovery <- compare_edge_recovery(
      original = true_model,
      simulated = model,
      threshold = 0.05
    )

    c(precision = recovery$precision,
      recall = recovery$recall,
      f1 = recovery$f1_score)
  })

  df <- do.call(rbind, recoveries)
  data.frame(
    n_sequences = n,
    precision_mean = mean(df[,"precision"]),
    recall_mean = mean(df[,"recall"]),
    f1_mean = mean(df[,"f1"])
  )
})

recovery_df <- do.call(rbind, recovery_results)
print(recovery_df)

# Typical power curve behavior
# - Steep improvement from 25-100
# - Diminishing returns after 200-300
# - Asymptotic approach to perfect recovery

# Compare different sequence lengths
length_grid <- generate_param_grid(
  param_ranges = list(
    n_sequences = c(50, 200),
    seq_length = c(10, 50),
    n_states = c(5, 5)
  ),
  n = 30,
  method = "lhs"
)

length_results <- run_grid_simulation(
  param_grid = length_grid,
  n_runs_per_setting = 30,
  model_type = "tna",
  seed = 42
)

# Compare different network sizes
states_grid <- generate_param_grid(
  param_ranges = list(
    n_sequences = c(50, 300),
    seq_length = c(25, 25),
    n_states = c(3, 10)
  ),
  n = 30,
  method = "lhs"
)

states_results <- run_grid_simulation(
  param_grid = states_grid,
  n_runs_per_setting = 30,
  model_type = "tna",
  seed = 42
)

# Test different network densities
# Use generate_probabilities with different parameters
sparse_net <- generate_probabilities(n_states = 6, edge_density = 0.3)
dense_net <- generate_probabilities(n_states = 6, edge_density = 0.8)

# 1. Define parameter space
full_grid <- generate_param_grid(
  param_ranges = list(
    n_sequences = c(25, 500),
    seq_length = c(15, 40),
    n_states = c(4, 8)
  ),
  n = 100,
  method = "lhs"
)

# 2. Run comprehensive simulation
full_results <- run_grid_simulation(
  param_grid = full_grid,
  n_runs_per_setting = 100,
  model_type = "tna",
  parallel = TRUE,
  cores = parallel::detectCores() - 1,
  seed = 42,
  verbose = TRUE
)

# 3. Analyze results
analysis <- summarize_grid_results(full_results)

# 4. Generate power tables
power_table <- summarize_simulation(
  full_results,
  by = c("n_sequences", "seq_length", "n_states"),
  metrics = c("mean", "sd", "ci_lower", "ci_upper")
)

# 5. Find recommended sample sizes
recommendations <- power_table[power_table$correlation_ci_lower >= 0.85, ]
min_by_condition <- aggregate(
  n_sequences ~ seq_length + n_states,
  data = recommendations,
  FUN = min
)
print(min_by_condition)

# Generate summary statistics
summary_stats <- summarize_simulation(
  power_results,
  by = "n_sequences",
  metrics = "all"
)

# Create report
cat("
Power Analysis Results
=====================

Study Parameters:
- Sequence length: 25
- Number of states: 5
- Replications per condition: 100
- Target correlation: 0.90

Results by Sample Size:
")
print(summary_stats)

cat("
Recommendations:
- Minimum sample size for r >= 0.90:", min_n, "sequences
- Recommended sample size (with margin):", ceiling(min_n * 1.2), "sequences
")