# Load models
MXL_wtp_Control_Not_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Control_Not_Pred")
MXL_wtp_Control_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Control_Pred")
MXL_wtp_Opt_Not_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Opt_Not_Pred")
MXL_wtp_Opt_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Opt_Pred")
MXL_wtp_Treated_Not_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Not_Pred")
MXL_wtp_Treated_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Pred")

# List of models
models <- list(
  "Treated Not Pred" = MXL_wtp_Treated_Not_Pred_model,
  "Control Not Pred" = MXL_wtp_Control_Not_Pred_model,
  "Control Pred" = MXL_wtp_Control_Pred_model,
  "Opt Not Pred" = MXL_wtp_Opt_Not_Pred_model,
  "Opt Pred" = MXL_wtp_Opt_Pred_model,
  "Treated Pred" = MXL_wtp_Treated_Pred_model
)

# Define the x-range for plotting (1000 values)
x_range <- seq(-20, 100, length.out = 1000)

# Define an empty data frame to store the results
plot_data <- data.frame()

# Loop through the models and calculate the density for each
for (i in seq_along(models)) {
  model_name <- names(models)[i]
  model <- models[[i]]
  
  # Extract mean and standard deviation for 'mu_natural' and 'sig_natural'
  coef <- model$estimate
  mean <- coef["mu_natural"]
  sd <- coef["sig_natural"]
  
  # Compute the normal distribution based on the extracted mean and sd
  y <- dnorm(x_range, mean, sd)
  
  # Create a data frame for this model and append it to the main data frame
  temp_data <- data.frame(x = x_range, density = y, model = model_name)
  plot_data <- bind_rows(plot_data, temp_data)
}

# Plot using ggplot2
ggplot(plot_data, aes(x = x, y = density, color = model)) +
  geom_line() + 
  labs(
    title = "Normal Distributions from Multiple MXL Models",
    x = "x",
    y = "Density",
    color = "Model"
  ) +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_color_brewer(palette = "Set1")

#### Z-test if distributions are different

z_value <- (MXL_wtp_Treated_Pred_model$estimate["mu_natural"] - MXL_wtp_Control_Pred_model$estimate["mu_natural"])/
  (sqrt((MXL_wtp_Treated_Pred_model$estimate["sig_natural"]/sqrt(MXL_wtp_Treated_Pred_model$nObs)) + 
     (MXL_wtp_Control_Pred_model$estimate["sig_natural"]/sqrt(MXL_wtp_Control_Pred_model$nObs))))


# Function to calculate z-value using your exact formula
calculate_z_value <- function(model_1, model_2) {
  # Calculate the z-value using your specific formula
  z_value <- (model_1$estimate["mu_natural"] - model_2$estimate["mu_natural"]) /
    (sqrt((model_1$estimate["sig_natural"] / sqrt(model_1$nObs/10)) + 
       (model_2$estimate["sig_natural"] / sqrt(model_2$nObs/10))))
  
  return(z_value)
}

# Example of using the function
z_value <- calculate_z_value(MXL_wtp_Treated_Pred_model, MXL_wtp_Control_Pred_model)

# Print the z-value
z_value