bootstrap_and_plot <- function(population, method_function, census=NA, n_replications = 1000) {
  
  # Capture the names of the dataset and method function
  dataset_name <- deparse(substitute(population))
  method_name <- deparse(substitute(method_function))
  
  data <- sample_population_small_area(population)
  
  # Bootstrap estimates
  start_time <- Sys.time()
  
  if (method_name %in% c("mrp", "lmp")) {
    replicated_means <- replicate(n_replications, method_function(data, census), simplify = FALSE)
  } else {
    replicated_means <- replicate(n_replications, method_function(data), simplify = FALSE)
  }
  
  duration <- round(difftime(Sys.time(), start_time, units = "secs"))
  print(paste("Execution time:", duration, "seconds."))
  
  # Combine all replicate means into one dataframe
  all_means <- bind_rows(replicated_means, .id = "replicate")
  
  # Calculate the overall mean and confidence intervals by small_area
  final_estimates <- all_means %>% 
    group_by(small_area) %>% 
    summarise(estimate_mean = mean(mean_attitude),
              lower_ci = quantile(mean_attitude, probs = 0.025),
              upper_ci = quantile(mean_attitude, probs = 0.975))
  
  # Calculate population averages for attitude by small_area in the original population
  population_averages <- population %>%
    group_by(small_area) %>%
    summarise(population_mean = mean(attitude))
  
  # Merge the population averages with the sampled estimates
  comparison_data <- final_estimates %>%
    inner_join(population_averages, by = "small_area") %>%
    mutate(error = estimate_mean - population_mean)
  
  # Calculate the proportion of small_areas where the confidence interval includes the population mean
  
  comparison_data <- comparison_data %>%
    mutate(in_ci = population_mean >= lower_ci & population_mean <= upper_ci)
  
  proportion_in_ci <- mean(comparison_data$in_ci)
  
  # Calculate MAE and Bias
  mae <- mean(abs(comparison_data$error))
  bias <- mean(comparison_data$error)
  
  # Calculate average variance across small_areas
  
  average_variance <- all_means %>% 
    group_by(small_area) %>% 
    summarise(variance = var(mean_attitude)) %>% 
    summarise(average_variance = mean(variance))
  
  # Plotting
  
  plot <- ggplot(comparison_data, aes(x = as.factor(small_area))) +
    geom_bar(aes(y = estimate_mean, fill = "Bootstrap Estimate"), stat = "identity", position = "dodge") +
    geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2, position = position_dodge(width = 0.9)) +
    geom_point(aes(y = population_mean, color = "Population Mean"), position = position_dodge(width = 0.9), size = 3) +
    theme_minimal() +
    labs(title = paste("Comparison of attitude Averages by small_area", "\nDataset:", dataset_name, "- Method:", method_name),
         subtitle = sprintf("MAE: %.2f | Bias: %.2f | Variance: %.2f | Prop in CI: %.2f", mae, bias, average_variance, proportion_in_ci),
         x = "small_area",
         y = "attitude Average") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    scale_fill_manual(name = "Estimate Type", values = c("Bootstrap Estimate" = "lightblue")) +
    scale_color_manual(name = "Estimate Type", values = c("Population Mean" = "darkred")) +
    ylim(0, 10)
  
  print(plot)
  
  # Define the file name
  file_name <- paste0("plot_", gsub(" ", "_", dataset_name), "_", gsub(" ", "_", method_name), ".png")
  
  # Save the plot
  ggsave(file_name, plot, width = 14, height = 8, units = "in")
}