diff --git a/.write-test-746fb7da-de5d-497a-baa2-d528cb118be7 b/.write-test-746fb7da-de5d-497a-baa2-d528cb118be7
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Scripts/logit/chr_vol_treat.R b/Scripts/logit/chr_vol_treat.R
index 5b494acb94e817a3895f194338df8a505a011609..caff63eda7317ef548129777474e943fb5f7fde4 100644
--- a/Scripts/logit/chr_vol_treat.R
+++ b/Scripts/logit/chr_vol_treat.R
@@ -109,78 +109,78 @@ testData <- labeled_data[-trainIndex, ]
 
 
 
-tuneGrid <- expand.grid(
-  nrounds = c(300),
-  max_depth = c(15),
-  eta = c(0.05),
-  gamma = c(0, 0.1),
-  colsample_bytree = c(0.6, 0.8),     # Fraction of features to be used for each tree
-  min_child_weight = c(1, 3, 5),      # Minimum sum of instance weight (hessian) needed in a child
-  subsample = c(0.7, 0.8)             # Fraction of samples to be used for each tree
-)
-
-
-model3 <- train(Choice_Treat ~ ., 
-                data = trainData, 
-                method = "xgbTree", 
-                tuneGrid = tuneGrid,
-                trControl = trainControl(method = "cv", number = 5))
-
-
-# Get variable importance
-varImp(model3)
-
-
-
-# Evaluate the model
-predictions <- predict(model3, newdata = testData)
-confusionMatrix(predictions, testData$Choice_Treat)
-labeled_predictions <- predict(model3, newdata = labeled_data)
-labeled_data$PredictedGroup <- labeled_predictions
-
-table(labeled_data$Choice_Treat, labeled_data$PredictedGroup)
-
-unlabeled_predictions <- predict(model3, newdata = unlabeled_data)
-labeled_data_id$PredictedGroup <- labeled_predictions
-data_prediction_labeled<-select(labeled_data_id, c("id", "PredictedGroup"))
-saveRDS(data_prediction_labeled, "Data/predictions_labeled.RDS")
-
-unlabeled_data$PredictedGroup <- unlabeled_predictions
-data_prediction<-select(unlabeled_data, c("id", "PredictedGroup"))
-saveRDS(data_prediction, "Data/predictions.RDS")
-
-print(model3$bestTune)
-
-
-table(data$predicted.classes ,data$Choice_Treat)
-
-calculate_metrics <- function(confusion_matrix) {
-  # Extract values from the confusion matrix
-  TN <- confusion_matrix[1, 1]
-  FP <- confusion_matrix[2, 1]
-  FN <- confusion_matrix[1, 2]
-  TP <- confusion_matrix[2, 2]
-  
-  # Calculate sensitivity
-  sensitivity <- round(TP / (TP + FN), 2)
-  
-  # Calculate specificity
-  specificity <- round(TN / (TN + FP), 2)
-  
-  # Return the results as a list
-  return(list(sensitivity = sensitivity, specificity = specificity))
-}
-
-# Calculate metrics
-metrics <- calculate_metrics(table(data$predicted.classes ,data$Choice_Treat))
-
-# Print results
-print(paste("Sensitivity:", metrics$sensitivity))
-print(paste("Specificity:", metrics$specificity))
-roc_obj <- roc(data$Choice_Treat, data$probabilities)
-plot(roc_obj, main = "ROC Curve", col = "blue", lwd = 2)
-auc_value <- auc(roc_obj)
-best_coords <- coords(roc_obj, "best", best.method="youden")
-
-
-cut_off <- best_coords$threshold
\ No newline at end of file
+# tuneGrid <- expand.grid(
+#   nrounds = c(300),
+#   max_depth = c(15),
+#   eta = c(0.05),
+#   gamma = c(0, 0.1),
+#   colsample_bytree = c(0.6, 0.8),     # Fraction of features to be used for each tree
+#   min_child_weight = c(1, 3, 5),      # Minimum sum of instance weight (hessian) needed in a child
+#   subsample = c(0.7, 0.8)             # Fraction of samples to be used for each tree
+# )
+# 
+# 
+# model3 <- train(Choice_Treat ~ ., 
+#                 data = trainData, 
+#                 method = "xgbTree", 
+#                 tuneGrid = tuneGrid,
+#                 trControl = trainControl(method = "cv", number = 5))
+# 
+# 
+# # Get variable importance
+# varImp(model3)
+# 
+# 
+# 
+# # Evaluate the model
+# predictions <- predict(model3, newdata = testData)
+# confusionMatrix(predictions, testData$Choice_Treat)
+# labeled_predictions <- predict(model3, newdata = labeled_data)
+# labeled_data$PredictedGroup <- labeled_predictions
+# 
+# table(labeled_data$Choice_Treat, labeled_data$PredictedGroup)
+# 
+# unlabeled_predictions <- predict(model3, newdata = unlabeled_data)
+# labeled_data_id$PredictedGroup <- labeled_predictions
+# data_prediction_labeled<-select(labeled_data_id, c("id", "PredictedGroup"))
+# saveRDS(data_prediction_labeled, "Data/predictions_labeled.RDS")
+# 
+# unlabeled_data$PredictedGroup <- unlabeled_predictions
+# data_prediction<-select(unlabeled_data, c("id", "PredictedGroup"))
+# saveRDS(data_prediction, "Data/predictions.RDS")
+# 
+# print(model3$bestTune)
+# 
+# 
+# table(data$predicted.classes ,data$Choice_Treat)
+# 
+# calculate_metrics <- function(confusion_matrix) {
+#   # Extract values from the confusion matrix
+#   TN <- confusion_matrix[1, 1]
+#   FP <- confusion_matrix[2, 1]
+#   FN <- confusion_matrix[1, 2]
+#   TP <- confusion_matrix[2, 2]
+#   
+#   # Calculate sensitivity
+#   sensitivity <- round(TP / (TP + FN), 2)
+#   
+#   # Calculate specificity
+#   specificity <- round(TN / (TN + FP), 2)
+#   
+#   # Return the results as a list
+#   return(list(sensitivity = sensitivity, specificity = specificity))
+# }
+# 
+# # Calculate metrics
+# metrics <- calculate_metrics(table(data$predicted.classes ,data$Choice_Treat))
+# 
+# # Print results
+# print(paste("Sensitivity:", metrics$sensitivity))
+# print(paste("Specificity:", metrics$specificity))
+# roc_obj <- roc(data$Choice_Treat, data$probabilities)
+# plot(roc_obj, main = "ROC Curve", col = "blue", lwd = 2)
+# auc_value <- auc(roc_obj)
+# best_coords <- coords(roc_obj, "best", best.method="youden")
+# 
+# 
+# cut_off <- best_coords$threshold
\ No newline at end of file
diff --git a/Scripts/logit/protesters.R b/Scripts/logit/protesters.R
index 722e38ff5be2444b3947a08e2cbf1fecce5c3f4c..bcbfbf80b8c037bf31455899cc97590dbcd7e80e 100644
--- a/Scripts/logit/protesters.R
+++ b/Scripts/logit/protesters.R
@@ -16,7 +16,7 @@ database_full <- database_full %>%
 
 data <- database_full %>%
   group_by(id) %>%
-  slice(1) %>%
+  dplyr::slice(1) %>%
   ungroup()
 
 data <- data %>% 
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Not_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Not_Pred.R
index 07ee9eb0b7d519a190714dd68a347e5f46f714ac..e31c769eddf9220d579df0d801a8f3a6c85c0895 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Not_Pred.R	
+++ b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Not_Pred.R	
@@ -142,7 +142,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 model = apollo_estimate(apollo_beta, apollo_fixed,
                         apollo_probabilities, apollo_inputs, 
                         estimate_settings=list(maxIterations=400,
-                                               estimationRoutine="bfgs",
+                                               estimationRoutine="bgw",
                                                hessianRoutine="analytic"))
 
 
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Not_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Not_Pred.R
index 484d1fd2c6bb631559313859713e7272a4bdd336..033604b4950ba7d6a7dee0f816440731cf6485e5 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Not_Pred.R	
+++ b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Not_Pred.R	
@@ -142,7 +142,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 model = apollo_estimate(apollo_beta, apollo_fixed,
                         apollo_probabilities, apollo_inputs, 
                         estimate_settings=list(maxIterations=400,
-                                               estimationRoutine="bfgs",
+                                               estimationRoutine="bgw",
                                                hessianRoutine="analytic"))
 
 
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Pred.R
index 0785b9d988270f61cc91e9529932869e2ae71a6e..92767fe966abe29e913a6568a83a0899005c0975 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Pred.R	
+++ b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Opt_Pred.R	
@@ -142,7 +142,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 model = apollo_estimate(apollo_beta, apollo_fixed,
                         apollo_probabilities, apollo_inputs, 
                         estimate_settings=list(maxIterations=400,
-                                               estimationRoutine="bfgs",
+                                               estimationRoutine="bgw",
                                                hessianRoutine="analytic"))
 
 
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Not_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Not_Pred.R
index 03667eaf757cf9ed466c2a5b1671e6944f8ec800..9ebe574cff14786ffc6aba841e96fe06a5e04105 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Not_Pred.R	
+++ b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Not_Pred.R	
@@ -142,7 +142,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 model = apollo_estimate(apollo_beta, apollo_fixed,
                         apollo_probabilities, apollo_inputs, 
                         estimate_settings=list(maxIterations=400,
-                                               estimationRoutine="bfgs",
+                                               estimationRoutine="bgw",
                                                hessianRoutine="analytic"))
 
 
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Pred.R
index 76b11f4fc54e4f86f0a6debe1532935b1565bfcd..dc5791d135f41ed5065319a92be2eaa98d76fc45 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Pred.R	
+++ b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Treatment_Pred.R	
@@ -142,7 +142,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 model = apollo_estimate(apollo_beta, apollo_fixed,
                                    apollo_probabilities, apollo_inputs, 
                                    estimate_settings=list(maxIterations=400,
-                                                          estimationRoutine="bfgs",
+                                                          estimationRoutine="bgw",
                                                           hessianRoutine="analytic"))
 
 
diff --git a/Scripts/ols/ols_consequentiality.R b/Scripts/ols/ols_consequentiality.R
index 7f89ba5e250506ce8c9f4e3bf0721278a111b689..71bb7d02450349a21120f741facb7398a86cb15d 100644
--- a/Scripts/ols/ols_consequentiality.R
+++ b/Scripts/ols/ols_consequentiality.R
@@ -1,35 +1,39 @@
-# Q28W3
-
-# Q29W3
-
-# A1 I belive, A5 I don't belive; A6 i don't know
-
-data <- data %>% mutate(Conseq_score = Conseq_UGS + Conseq_Money)
-
-
-conseq_model_A <- lm(Conseq_score ~ as.factor(Treatment_A), data) 
-summary(conseq_model_A)
-
-conseq_model_control_A <- lm(Conseq_score ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(conseq_model_control_A )
-
-
-conseq_model_B <- lm(Conseq_score ~ as.factor(Treatment_B), data) 
-summary(conseq_model_B)
-
-conseq_model_control_B <- lm(Conseq_score ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(conseq_model_control_B)
-
-
-conseq_model_C <- lm(Conseq_score ~ as.factor(Treatment_C), data) 
-summary(conseq_model_C)
-
-conseq_model_control_C <- lm(Conseq_score ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(conseq_model_control_C)
-
-
-conseq_model_D <- lm(Conseq_score ~ as.factor(Treatment_D), data) 
-summary(conseq_model_D)
-
-conseq_model_control_D <- lm(Conseq_score ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(conseq_model_control_D)
+# Q28W3
+
+# Q29W3
+
+# A1 I belive, A5 I don't belive; A6 i don't know
+
+data <- data %>% mutate(Conseq_score = Conseq_UGS + Conseq_Money)
+
+
+conseq_model_A <- lm(Conseq_score ~ as.factor(Treatment_A), data) 
+summary(conseq_model_A)
+
+conseq_model_control_A <- lm(Conseq_score ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(conseq_model_control_A )
+
+
+conseq_model_B <- lm(Conseq_score ~ as.factor(Treatment_B), data) 
+summary(conseq_model_B)
+
+conseq_model_control_B <- lm(Conseq_score ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(conseq_model_control_B)
+
+
+conseq_model_C <- lm(Conseq_score ~ as.factor(Treatment_C), data) 
+summary(conseq_model_C)
+
+conseq_model_control_C <- lm(Conseq_score ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(conseq_model_control_C)
+
+
+conseq_model_D <- lm(Conseq_score ~ as.factor(Treatment_D), data) 
+summary(conseq_model_D)
+
+conseq_model_control_D <- lm(Conseq_score ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(conseq_model_control_D)
+
+
+conseq_model_control_G <- lm(Conseq_score ~ as.factor(Matching_Group) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(conseq_model_control_G)
diff --git a/Scripts/ols/ols_opt_out.R b/Scripts/ols/ols_opt_out.R
index 74784e9ab3268b04ac2c67a5b830bfb8e9e704dd..64e5dffd972444583f08107d1ce434453fc32e07 100644
--- a/Scripts/ols/ols_opt_out.R
+++ b/Scripts/ols/ols_opt_out.R
@@ -1,87 +1,87 @@
-data <- database_full %>%
-  group_by(id) %>%
-  slice(1) %>%
-  ungroup()
-
-data$Treatment_C <- as.factor(data$Treatment_C)
-
-data$Treatment_C <- relevel(data$Treatment_C, ref = "No Treatment 3")
-data$Treatment_D <- as.factor(data$Treatment_D)
-
-data$Treatment_D <- relevel(data$Treatment_D, ref = "No Treatment 3")
-ols_opt_out_A<- lm( count_choosen_3 ~ as.factor(Treatment_A) ,data)
-summary(ols_opt_out_A)
-ols_opt_out_control_A<- lm( count_choosen_3 ~ as.factor(Treatment_A) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
-summary(ols_opt_out_control_A)
-
-ols_opt_out_B<- lm( count_choosen_3 ~ as.factor(Treatment_B) ,data)
-summary(ols_opt_out_B)
-ols_opt_out_control_B<- lm( count_choosen_3 ~ as.factor(Treatment_B) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
-summary(ols_opt_out_control_B)
-ols_opt_out_C<- lm( count_choosen_3 ~ as.factor(Treatment_C) ,data)
-summary(ols_opt_out_C)
-ols_opt_out_control_C<- lm( count_choosen_3 ~ as.factor(Treatment_C) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
-summary(ols_opt_out_control_C)
-ols_opt_out_D<- lm( count_choosen_3 ~ as.factor(Treatment_D) ,data)
-summary(ols_opt_out_D)
-ols_opt_out_control_D<- lm( count_choosen_3 ~ as.factor(Treatment_D) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
-summary(ols_opt_out_control_D)
-
-
-# Obtain predicted values
-predictions <- predict(ols_opt_out_control_C, data, se.fit = TRUE)
-
-# Create a data frame with predictions and standard errors
-predictions_df <- data.frame(
-  Treatment_C = data$Treatment_C,
-  Predicted_Value = predictions$fit,
-  Standard_Error = predictions$se.fit
-)
-
-# Calculate means and standard errors
-means <- aggregate(Predicted_Value ~ Treatment_C, predictions_df, mean)
-means$Standard_Error <- aggregate(Standard_Error ~ Treatment_C, predictions_df, mean)$Standard_Error
-
-# Plot mean and standard errors
-library(ggplot2)
-
-ggplot(means, aes(x = as.factor(Treatment_C), y = Predicted_Value, fill = as.factor(Treatment_C))) +
-  geom_bar(stat = "identity", position = "dodge") +
-  geom_errorbar(aes(ymin = Predicted_Value - 1.64*Standard_Error, ymax = Predicted_Value + 1.64*Standard_Error),
-                position = position_dodge(width = 0.9), width = 0.25) +
-  labs(fill = "Treatment") +
-  xlab("Treatment") +
-  ylab("Mean Predicted Value") +
-  theme_minimal()
-# 
-# # Create an HTML results table with customized names and stars
-# results_table_1 <- stargazer(
-#   ols_opt_out_1, ols_opt_out_control_1, 
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "No. Opt-out choices",
-#   covariate.labels = c("Video 1", "Video 2", "No Info choosen", "Text 1", "Text 2",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_1, "Estimation_results/ols/opt_out_1.html")
-# 
-# # Create an HTML results table with customized names and stars
-# results_table_2 <- stargazer(
-#   ols_opt_out_2, ols_opt_out_control_2, 
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "No. Opt-out choices",
-#   covariate.labels = c("Video", "No Info choosen", "Text",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_2, "Estimation_results/ols/opt_out_2.html")
+data <- database_full %>%
+  group_by(id) %>%
+  dplyr::slice(1) %>%
+  ungroup()
+
+data$Treatment_C <- as.factor(data$Treatment_C)
+
+data$Treatment_C <- relevel(data$Treatment_C, ref = "No Treatment 3")
+data$Treatment_D <- as.factor(data$Treatment_D)
+
+data$Treatment_D <- relevel(data$Treatment_D, ref = "No Treatment 3")
+ols_opt_out_A<- lm( count_choosen_3 ~ as.factor(Treatment_A) ,data)
+summary(ols_opt_out_A)
+ols_opt_out_control_A<- lm( count_choosen_3 ~ as.factor(Treatment_A) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
+summary(ols_opt_out_control_A)
+
+ols_opt_out_B<- lm( count_choosen_3 ~ as.factor(Treatment_B) ,data)
+summary(ols_opt_out_B)
+ols_opt_out_control_B<- lm( count_choosen_3 ~ as.factor(Treatment_B) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
+summary(ols_opt_out_control_B)
+ols_opt_out_C<- lm( count_choosen_3 ~ as.factor(Treatment_C) ,data)
+summary(ols_opt_out_C)
+ols_opt_out_control_C<- lm( count_choosen_3 ~ as.factor(Treatment_C) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
+summary(ols_opt_out_control_C)
+ols_opt_out_D<- lm( count_choosen_3 ~ as.factor(Treatment_D) ,data)
+summary(ols_opt_out_D)
+ols_opt_out_control_D<- lm( count_choosen_3 ~ as.factor(Treatment_D) + Z_Mean_NR + QFIncome + as.factor(Gender)+Age_mean+Uni_degree,data)
+summary(ols_opt_out_control_D)
+
+
+# Obtain predicted values
+predictions <- predict(ols_opt_out_control_C, data, se.fit = TRUE)
+
+# Create a data frame with predictions and standard errors
+predictions_df <- data.frame(
+  Treatment_C = data$Treatment_C,
+  Predicted_Value = predictions$fit,
+  Standard_Error = predictions$se.fit
+)
+
+# Calculate means and standard errors
+means <- aggregate(Predicted_Value ~ Treatment_C, predictions_df, mean)
+means$Standard_Error <- aggregate(Standard_Error ~ Treatment_C, predictions_df, mean)$Standard_Error
+
+# Plot mean and standard errors
+library(ggplot2)
+
+ggplot(means, aes(x = as.factor(Treatment_C), y = Predicted_Value, fill = as.factor(Treatment_C))) +
+  geom_bar(stat = "identity", position = "dodge") +
+  geom_errorbar(aes(ymin = Predicted_Value - 1.64*Standard_Error, ymax = Predicted_Value + 1.64*Standard_Error),
+                position = position_dodge(width = 0.9), width = 0.25) +
+  labs(fill = "Treatment") +
+  xlab("Treatment") +
+  ylab("Mean Predicted Value") +
+  theme_minimal()
+# 
+# # Create an HTML results table with customized names and stars
+# results_table_1 <- stargazer(
+#   ols_opt_out_1, ols_opt_out_control_1, 
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "No. Opt-out choices",
+#   covariate.labels = c("Video 1", "Video 2", "No Info choosen", "Text 1", "Text 2",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_1, "Estimation_results/ols/opt_out_1.html")
+# 
+# # Create an HTML results table with customized names and stars
+# results_table_2 <- stargazer(
+#   ols_opt_out_2, ols_opt_out_control_2, 
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "No. Opt-out choices",
+#   covariate.labels = c("Video", "No Info choosen", "Text",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_2, "Estimation_results/ols/opt_out_2.html")
diff --git a/Scripts/ols/ols_quiz.R b/Scripts/ols/ols_quiz.R
index 489286f049b6757b9215747ef045532ae70a0251..176f0e517b808c6e5c8edf6078003e7c570a502d 100644
--- a/Scripts/ols/ols_quiz.R
+++ b/Scripts/ols/ols_quiz.R
@@ -1,70 +1,77 @@
-
-quiz_data$Treatment_C <- as.factor(quiz_data$Treatment_C)
-quiz_data$Treatment_D <- as.factor(quiz_data$Treatment_D)
-
-quiz_data$Treatment_C <- relevel(quiz_data$Treatment_C, ref = "No Treatment 3")
-quiz_data$Treatment_D <- relevel(quiz_data$Treatment_D, ref = "No Treatment 3")
-
-ols_percentage_correct_A<- lm( percentage_correct ~ as.factor(Treatment_A) ,quiz_data)
-summary(ols_percentage_correct_A)
-ols_percentage_correct_control_A<- lm( percentage_correct ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
-summary(ols_percentage_correct_control_A)
-
-ols_percentage_correct_B<- lm( percentage_correct ~ as.factor(Treatment_B) ,quiz_data)
-summary(ols_percentage_correct_B)
-ols_percentage_correct_control_B<- lm( percentage_correct ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
-summary(ols_percentage_correct_control_B)
-
-ols_percentage_correct_C<- lm( percentage_correct ~ as.factor(Treatment_C) ,quiz_data)
-summary(ols_percentage_correct_C)
-ols_percentage_correct_control_C<- lm( percentage_correct ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
-summary(ols_percentage_correct_control_C)
-
-vif(ols_percentage_correct_control_C)
-
-
-ols_percentage_correct_D<- lm( percentage_correct ~ as.factor(Treatment_D) ,quiz_data)
-summary(ols_percentage_correct_D)
-ols_percentage_correct_control_D<- lm( percentage_correct ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
-summary(ols_percentage_correct_control_D)
-# 
-# # Create an HTML results table with customized names and stars
-# results_table_3 <- stargazer(
-#   ols_quiz_1a, ols_quiz_control_1a, ols_quiz_1b, ols_quiz_control_1b,  
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "No. correct quiz questions",
-#   covariate.labels = c("Video 1", "Video 2", "No Info choosen", "Text 1", "Text 2",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_3, "Estimation_results/ols/quiz_1.html")
-# 
-# # Create an HTML results table with customized names and stars
-# results_table_4 <- stargazer(
-#   ols_quiz_2a, ols_quiz_control_2a, ols_quiz_2b, ols_quiz_control_2b,  
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "No. correct quiz questions",
-#   covariate.labels = c("Video", "No Info choosen", "Text",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#     star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_4, "Estimation_results/ols/quiz_2.html")
-# 
-# 
-# # Poisson regression
-# poisson_quiz_2b <- glm(
-#   number_correct ~ Dummy_Video + Dummy_no_info + Dummy_nv,
-#   quiz_data = quiz_data,
-#   family = poisson
-# )
-# summary(poisson_quiz_2b)
+matching_group <- data %>% dplyr::select(id, Matching_Group)
+
+quiz_data <- left_join(quiz_data, matching_group, by="id")
+
+
+quiz_data$Treatment_C <- as.factor(quiz_data$Treatment_C)
+quiz_data$Treatment_D <- as.factor(quiz_data$Treatment_D)
+
+quiz_data$Treatment_C <- relevel(quiz_data$Treatment_C, ref = "No Treatment 3")
+quiz_data$Treatment_D <- relevel(quiz_data$Treatment_D, ref = "No Treatment 3")
+
+ols_percentage_correct_A<- lm( percentage_correct ~ as.factor(Treatment_A) ,quiz_data)
+summary(ols_percentage_correct_A)
+ols_percentage_correct_control_A<- lm( percentage_correct ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
+summary(ols_percentage_correct_control_A)
+
+ols_percentage_correct_B<- lm( percentage_correct ~ as.factor(Treatment_B) ,quiz_data)
+summary(ols_percentage_correct_B)
+ols_percentage_correct_control_B<- lm( percentage_correct ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
+summary(ols_percentage_correct_control_B)
+
+ols_percentage_correct_C<- lm( percentage_correct ~ as.factor(Treatment_C) ,quiz_data)
+summary(ols_percentage_correct_C)
+ols_percentage_correct_control_C<- lm( percentage_correct ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
+summary(ols_percentage_correct_control_C)
+
+vif(ols_percentage_correct_control_C)
+
+
+ols_percentage_correct_D<- lm( percentage_correct ~ as.factor(Treatment_D) ,quiz_data)
+summary(ols_percentage_correct_D)
+ols_percentage_correct_control_D<- lm( percentage_correct ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
+summary(ols_percentage_correct_control_D)
+
+ols_percentage_correct_control_G<- lm( percentage_correct ~ as.factor(Matching_Group) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,quiz_data)
+summary(ols_percentage_correct_control_G)
+# 
+# # Create an HTML results table with customized names and stars
+# results_table_3 <- stargazer(
+#   ols_quiz_1a, ols_quiz_control_1a, ols_quiz_1b, ols_quiz_control_1b,  
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "No. correct quiz questions",
+#   covariate.labels = c("Video 1", "Video 2", "No Info choosen", "Text 1", "Text 2",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_3, "Estimation_results/ols/quiz_1.html")
+# 
+# # Create an HTML results table with customized names and stars
+# results_table_4 <- stargazer(
+#   ols_quiz_2a, ols_quiz_control_2a, ols_quiz_2b, ols_quiz_control_2b,  
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "No. correct quiz questions",
+#   covariate.labels = c("Video", "No Info choosen", "Text",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#     star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_4, "Estimation_results/ols/quiz_2.html")
+# 
+# 
+# # Poisson regression
+# poisson_quiz_2b <- glm(
+#   number_correct ~ Dummy_Video + Dummy_no_info + Dummy_nv,
+#   quiz_data = quiz_data,
+#   family = poisson
+# )
+# summary(poisson_quiz_2b)
diff --git a/Scripts/ols/ols_time_spent.R b/Scripts/ols/ols_time_spent.R
index b417979414bcea8a8105ba80531ee72a85faa6dc..fe6de07bf090047c18e4c6cb41b01a77c27630e8 100644
--- a/Scripts/ols/ols_time_spent.R
+++ b/Scripts/ols/ols_time_spent.R
@@ -1,79 +1,110 @@
-library(car)
-
-data <- database_full %>%
-  group_by(id) %>%
-  dplyr::slice(1) %>%
-  ungroup()
-
-
-data$Treatment_C <- as.factor(data$Treatment_C)
-data$Treatment_D <- as.factor(data$Treatment_D)
-data$Treatment_C <- relevel(data$Treatment_C, ref = "No Treatment 3")
-data$Treatment_D <- relevel(data$Treatment_D, ref = "No Treatment 3")
-ols_time_spent_A<- lm( interviewtime_net_clean ~ as.factor(Treatment_A) ,data)
-summary(ols_time_spent_A)
-ols_time_spent_control_A<- lm( interviewtime_net_clean ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(ols_time_spent_control_A)
-
-ols_time_spent_B<- lm( interviewtime_net_clean ~ as.factor(Treatment_B) ,data)
-summary(ols_time_spent_B)
-ols_time_spent_control_B<- lm( interviewtime_net_clean ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
-summary(ols_time_spent_control_B)
-ols_time_spent_C<- lm( interviewtime_net_clean ~ as.factor(Treatment_C) ,data)
-summary(ols_time_spent_C)
-ols_time_spent_control_C<- lm( interviewtime_net_clean ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(ols_time_spent_control_C)
-ols_time_spent_D<- lm( interviewtime_net_clean ~ as.factor(Treatment_D) ,data)
-summary(ols_time_spent_D)
-ols_time_spent_control_D<- lm( interviewtime_net_clean ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
-summary(ols_time_spent_control_D)
-
-ols_time_cc_A<- lm( CC_time_mean_clean ~ as.factor(Treatment_A) ,data)
-summary(ols_time_cc_A)
-ols_time_cc_control_A<- lm( CC_time_mean_clean ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
-summary(ols_time_cc_control_A)
-
-ols_time_cc_B<- lm( CC_time_mean_clean ~ as.factor(Treatment_B) ,data)
-summary(ols_time_cc_B)
-ols_time_cc_control_B<- lm( CC_time_mean_clean ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome + Uni_degree,data)
-summary(ols_time_cc_control_B)
-ols_time_cc_C<- lm( CC_time_mean_clean ~ as.factor(Treatment_C) ,data)
-summary(ols_time_cc_C)
-ols_time_cc_control_C<- lm( CC_time_mean_clean ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
-summary(ols_time_cc_control_C)
-ols_time_cc_D<- lm( CC_time_mean_clean ~ as.factor(Treatment_D) ,data)
-summary(ols_time_cc_D)
-ols_time_cc_control_D<- lm( CC_time_mean_clean ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
-summary(ols_time_cc_control_D)
-# # Create an HTML results table with customized names and stars
-# results_table_5 <- stargazer(
-#   ols_tme_spent_1, ols_tme_spent_control_1, 
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "Interview Time (without Treatments)",
-#   covariate.labels = c("Treated_A", "Treated_B","No Info choosen", "Video 1", "Video 2",  "Text 1", "Text 2",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_5, "Estimation_results/ols/time_spent_1.html")
-# 
-# # Create an HTML results table with customized names and stars
-# results_table_6 <- stargazer(
-#   ols_tme_spent_2, ols_tme_spent_control_2, 
-#   align = TRUE,
-#   type = "html",
-#   dep.var.labels = "Interview Time (without Treatments)",
-#   covariate.labels = c("Video", "No Info choosen", "Text",
-#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
-#                        "Hochschulabschluss") ,
-#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
-#   star.char = c("*", "**", "***")  # Custom significance stars
-# )
-# 
-# # Save the HTML table to a file
-# writeLines(results_table_6, "Estimation_results/ols/time_spent_2.html")
-
+library(car)
+
+
+data_predictions1 <- readRDS("Data/predictions.RDS")
+data_predictions2 <- readRDS("Data/predictions_labeled.RDS")
+
+data_predictions <- bind_rows(data_predictions1, data_predictions2)
+
+database <- left_join(database_full, data_predictions, by="id")
+
+
+
+database <- database %>% 
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0)) %>% 
+  mutate(Matching_Group = case_when(Dummy_Treated == 1 & PredictedGroup == 1 ~"Treat Pred",
+                                    Dummy_Treated == 1 & PredictedGroup == 0 ~"Treat Not Pred",
+                                    Treatment_A == "Vol_Treated" & PredictedGroup == 1 ~ "Opt Pred",
+                                    Treatment_A == "Vol_Treated" & PredictedGroup == 0 ~ "Opt Not Pred",
+                                    Treatment_new == 6 & PredictedGroup == 0 ~ "Control Pred",
+                                    TRUE ~ "Control Not Pred"))
+
+data <- database %>%
+  group_by(id) %>%
+  dplyr::slice(1) %>%
+  ungroup()
+
+
+data$Treatment_C <- as.factor(data$Treatment_C)
+data$Treatment_D <- as.factor(data$Treatment_D)
+data$Treatment_C <- relevel(data$Treatment_C, ref = "No Treatment 3")
+data$Treatment_D <- relevel(data$Treatment_D, ref = "No Treatment 3")
+ols_time_spent_A<- lm( interviewtime_net_clean ~ as.factor(Treatment_A) ,data)
+summary(ols_time_spent_A)
+ols_time_spent_control_A<- lm( interviewtime_net_clean ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(ols_time_spent_control_A)
+
+ols_time_spent_B<- lm( interviewtime_net_clean ~ as.factor(Treatment_B) ,data)
+summary(ols_time_spent_B)
+ols_time_spent_control_B<- lm( interviewtime_net_clean ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
+summary(ols_time_spent_control_B)
+ols_time_spent_C<- lm( interviewtime_net_clean ~ as.factor(Treatment_C) ,data)
+summary(ols_time_spent_C)
+ols_time_spent_control_C<- lm( interviewtime_net_clean ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(ols_time_spent_control_C)
+ols_time_spent_D<- lm( interviewtime_net_clean ~ as.factor(Treatment_D) ,data)
+summary(ols_time_spent_D)
+ols_time_spent_control_D<- lm( interviewtime_net_clean ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(ols_time_spent_control_D)
+
+ols_time_spent_control_G<- lm( interviewtime_net_clean ~ as.factor(Matching_Group) + Z_Mean_NR + as.factor(Gender)+Age_mean+ QFIncome + Uni_degree,data)
+summary(ols_time_spent_control_G)
+
+
+
+ols_time_cc_A<- lm( CC_time_mean_clean ~ as.factor(Treatment_A) ,data)
+summary(ols_time_cc_A)
+ols_time_cc_control_A<- lm( CC_time_mean_clean ~ as.factor(Treatment_A) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
+summary(ols_time_cc_control_A)
+
+ols_time_cc_B<- lm( CC_time_mean_clean ~ as.factor(Treatment_B) ,data)
+summary(ols_time_cc_B)
+ols_time_cc_control_B<- lm( CC_time_mean_clean ~ as.factor(Treatment_B) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome + Uni_degree,data)
+summary(ols_time_cc_control_B)
+ols_time_cc_C<- lm( CC_time_mean_clean ~ as.factor(Treatment_C) ,data)
+summary(ols_time_cc_C)
+ols_time_cc_control_C<- lm( CC_time_mean_clean ~ as.factor(Treatment_C) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
+summary(ols_time_cc_control_C)
+ols_time_cc_D<- lm( CC_time_mean_clean ~ as.factor(Treatment_D) ,data)
+summary(ols_time_cc_D)
+ols_time_cc_control_D<- lm( CC_time_mean_clean ~ as.factor(Treatment_D) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
+summary(ols_time_cc_control_D)
+
+
+ols_time_cc_control_G<- lm( CC_time_mean_clean ~ as.factor(Matching_Group) + Z_Mean_NR + as.factor(Gender)+Age_mean + QFIncome +Uni_degree,data)
+summary(ols_time_cc_control_G)
+# # Create an HTML results table with customized names and stars
+# results_table_5 <- stargazer(
+#   ols_tme_spent_1, ols_tme_spent_control_1, 
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "Interview Time (without Treatments)",
+#   covariate.labels = c("Treated_A", "Treated_B","No Info choosen", "Video 1", "Video 2",  "Text 1", "Text 2",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_5, "Estimation_results/ols/time_spent_1.html")
+# 
+# # Create an HTML results table with customized names and stars
+# results_table_6 <- stargazer(
+#   ols_tme_spent_2, ols_tme_spent_control_2, 
+#   align = TRUE,
+#   type = "html",
+#   dep.var.labels = "Interview Time (without Treatments)",
+#   covariate.labels = c("Video", "No Info choosen", "Text",
+#                        "Female", "Divers", "Age", "mittlere Reife", "Abitur", "Berufsausbildung",
+#                        "Hochschulabschluss") ,
+#   star.cutoffs = c(0.1, 0.05, 0.01),  # Custom significance levels
+#   star.char = c("*", "**", "***")  # Custom significance stars
+# )
+# 
+# # Save the HTML table to a file
+# writeLines(results_table_6, "Estimation_results/ols/time_spent_2.html")
+
diff --git a/Scripts/visualize_distr_spli.R b/Scripts/visualize_distr_spli.R
index e39949caa4bacb1d93ee429b8774cd4347094536..09818a06f89996c90f08c3097fea5fdbddd7350c 100644
--- a/Scripts/visualize_distr_spli.R
+++ b/Scripts/visualize_distr_spli.R
@@ -1,10 +1,10 @@
 # Load models
-MXL_wtp_Control_Not_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Control_Not_Pred_model.rds")
-MXL_wtp_Control_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Control_Pred_model.rds")
-MXL_wtp_Opt_Not_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Opt_Not_Pred_model.rds")
-MXL_wtp_Opt_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Opt_Pred_model.rds")
-MXL_wtp_Treated_Not_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Not_Pred_model.rds")
-MXL_wtp_Treated_Pred_model <- readRDS("C:/nextcloud/Hot_Topic_Cool_Choices/Estimation_results/mxl/prediction/Split_samples/MXL_wtp_Treated_Pred_model.rds")
+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(
@@ -17,39 +17,62 @@ models <- list(
 )
 
 # Define the x-range for plotting (1000 values)
-x_range <- seq(-10, 100, length.out = 1000)
+x_range <- seq(-20, 100, length.out = 1000)
 
-# Define colors and line types for each model
-colors <- c("blue", "red", "green", "purple", "orange", "brown")
-line_types <- 1:6
+# Define an empty data frame to store the results
+plot_data <- data.frame()
 
-# Initialize the plot with the first model
-model <- models[["Treated Not Pred"]]
-coef <- summary(model)$estimate
-mean <- coef["mu_natural", "Estimate"]
-sd <- coef["sig_natural", "Estimate"]
-y <- dnorm(x_range, mean, sd)
-
-plot(x_range, y, type = "l", col = colors[1], lwd = 2, lty = line_types[1],
-     ylim = c(0, max(y)), xlab = "x", ylab = "Density", 
-     main = "Normal Distributions from Multiple MXL Models")
-
-# Loop through the remaining models and add their distributions to the plot
-for (i in 2:length(models)) {
+# 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 <- summary(model)$estimate
-  mean <- coef["mu_natural", "Estimate"]
-  sd <- coef["sig_natural", "Estimate"]
+  coef <- model$estimate
+  mean <- coef["mu_natural"]
+  sd <- coef["sig_natural"]
   
-  # Compute the normal distribution based on extracted mean and sd
+  # Compute the normal distribution based on the extracted mean and sd
   y <- dnorm(x_range, mean, sd)
   
-  # Add the distribution to the existing plot
-  lines(x_range, y, col = colors[i], lwd = 2, lty = line_types[i])
+  # 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)
 }
 
-# Add a legend
-legend("topright", inset = c(-0.25, 0), legend = names(models), col = colors, 
-       lty = line_types, lwd = 2, xpd = TRUE)
+# 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
+