PSM start

Scripts/MAKE_FILE.R

@@ -107,3 +107,14 @@ source("Scripts/interaction_plots_presi.R")
 # mxl_wtp_case_b_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp Case B prot")
 # mxl_wtp_case_c_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp_Case_C prot")
 
+
+#### Implement Propensity Score Matching 
+
+
+####### Show distributions of split samples and compare 
+
+source("Scripts/visualize_distr_spli.R")
+
+
+
+

Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Pred.R

@@ -143,7 +143,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"))
 
 

Scripts/ols/ols_consequentiality.R

@@ -35,5 +35,5 @@ conseq_model_control_D <- lm(Conseq_score ~ as.factor(Treatment_D) + Z_Mean_NR +
 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)
+# 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)

Scripts/propensity_score_matching.R

+library(MatchIt)
+library(cobalt)
+
+data_psm_control_treat <- data %>% select_if(~ !all(is.na(.))) %>%
+  select(Treatment_A, Z_Mean_NR, Age, Uni_degree, Kids_Dummy, Gender_female, Naturalness_SQ) %>% 
+  filter(Treatment_A != "Vol_Treated") %>% mutate(Dummy_Treat = case_when(Treatment_A == "Treated" ~1, TRUE~0))
+
+
+psm_model <- matchit(Dummy_Treat ~ Z_Mean_NR + Age + Uni_degree + Kids_Dummy + Gender_female + Naturalness_SQ,
+                     data=data_psm_control_treat, method="full",
+                     distance="glm", ratio=1)
+
+
+psm_summary <- summary(psm_model)
+print(psm_summary)
+
+
+matched_data <- match.data(psm_model)
+
+love.plot(psm_model, stats = "mean.diffs", abs = FALSE, drop.distance = FALSE, 
+          thresholds = c(m = 0.1))
+
+
+bal.plot(psm_model, var.name = "distance", which = "both", type = "density")
+

Scripts/visualize_distr_spli.R

@@ -6,6 +6,11 @@ MXL_wtp_Opt_Pred_model <- apollo_loadModel("Estimation_results/mxl/prediction/Sp
 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")
 
+
+mxl_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Treated A")
+mxl_vol_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Treated A")
+mxl_not_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Not_Treated A")
+
 # List of models
 models <- list(
   "Treated Not Pred" = MXL_wtp_Treated_Not_Pred_model,
@@ -13,19 +18,23 @@ models <- list(
   "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
+  "Treated Pred" = MXL_wtp_Treated_Pred_model,
+  "Control Split Sample" = mxl_not_tr,
+  "Opt Split Sample" = mxl_vol_tr,
+  "Treated Split Sample" = mxl_tr
 )
 
 group_color_mapping <- c(
   "Control" = "grey",
   "Treated" = "red",
-  "Optional" = "yellow"
+  "Optional" = "orange"
 )
 
 # Define line type mapping for prediction status
 prediction_line_mapping <- c(
   "Not Predicted" = "dashed",
-  "Predicted" = "solid"
+  "Predicted" = "solid",
+  "Split Sample" = "dotted"
 )
 
 # Create a mapping for the groups
@@ -35,17 +44,23 @@ group_mapping <- c(
   "Treated Not Pred" = "Treated",
   "Treated Pred" = "Treated",
   "Opt Not Pred" = "Optional",
-  "Opt Pred" = "Optional"
+  "Opt Pred" = "Optional",
+  "Control Split Sample" = "Control",
+  "Opt Split Sample" = "Optional",
+  "Treated Split Sample" = "Treated"
 )
 
 # Create a mapping for prediction status
 prediction_mapping <- c(
   "Control Not Pred" = "Not Predicted",
   "Control Pred" = "Predicted",
+  "Control Split Sample" = "Split Sample",
   "Treated Not Pred" = "Not Predicted",
   "Treated Pred" = "Predicted",
+  "Treated Split Sample" = "Split Sample",
   "Opt Not Pred" = "Not Predicted",
-  "Opt Pred" = "Predicted"
+  "Opt Pred" = "Predicted",
+  "Opt Split Sample" = "Split Sample"
 )
 
 # Define the x-range for plotting (1000 values)
@@ -117,4 +132,87 @@ z_value <- calculate_z_value(MXL_wtp_Treated_Pred_model, MXL_wtp_Control_Pred_mo
 
 # Print the z-value
 z_value
-   
+
+
+#### Plot underneath
+
+ggplot(plot_data, aes(x = x, y = density, color = Group, linetype = Prediction)) +
+  geom_line(size = 1) + 
+  scale_color_manual(values = group_color_mapping) +  # Set custom colors for groups
+  scale_linetype_manual(values = prediction_line_mapping) +  # Set custom linetypes for prediction status
+  facet_grid(Prediction ~ ., scales = "free_y") +  # Stacks "Not Predicted" on top of "Predicted"
+  labs(
+    title = "Normal Distributions from Multiple MXL Models",
+    x = "WTP Naturalness (€/month)",
+    y = "Density",
+    color = "Group",
+    linetype = "Prediction"
+  ) +
+  theme_minimal() +
+  theme(
+    legend.position = "right",
+    strip.text.y = element_text(size = 12, face = "bold"),  # Adjust facet labels
+    panel.spacing = unit(1, "lines")  # Adds space between the two facets
+  ) +
+  xlim(-25, 75)
+
+
+#### Compare means
+# Create an empty data frame
+# Function to compute CI dynamically
+calculate_ci <- function(mean, se, confidence = 0.90) {
+  z_value <- qnorm((1 + confidence) / 2)  # Compute Z-score for the given confidence level
+  lower_ci <- mean - z_value * se
+  upper_ci <- mean + z_value * se
+  return(c(lower_ci, upper_ci))
+}
+
+ci_data <- data.frame()
+
+# Loop through each model and extract mean and confidence interval
+for (i in seq_along(models)) {
+  model_name <- names(models)[i]
+  model <- models[[i]]
+  
+  # Extract mean estimate and robust standard error
+  mean_value <- model$estimate["mu_natural"]
+  se_robust <- model$robse["mu_natural"]
+  
+  # Compute 90% confidence intervals dynamically
+  ci_values <- calculate_ci(mean_value, se_robust, confidence = 0.90)
+  
+  # Append to data frame
+  temp_data <- data.frame(
+    Model = model_name,
+    Mean = mean_value,
+    Lower_CI = ci_values[1],
+    Upper_CI = ci_values[2],
+    Group = group_mapping[[model_name]], 
+    Prediction = prediction_mapping[[model_name]]
+  )
+  
+  ci_data <- bind_rows(ci_data, temp_data)
+}
+
+# Ensure correct order of Groups (Treated > Optional > Control)
+ci_data$Group <- factor(ci_data$Group, levels = c("Treated", "Optional", "Control"))
+ci_data$Model <- factor(ci_data$Model, levels = unique(ci_data$Model))
+
+ggplot(ci_data, aes(x = Model, y = Mean, color = Group)) +
+  geom_point(size = 3, position = position_dodge(width = 0.5)) +  # Mean points
+  geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2, position = position_dodge(width = 0.5)) +  # CI bars
+  scale_color_manual(values = group_color_mapping) +  # Custom colors for groups
+  facet_grid(Prediction ~ ., scales = "free_y") +  # Stack "Not Predicted" above "Predicted"
+  coord_flip() +  # Flip axes for better readability
+  labs(
+    title = "Mean WTP Naturalness Estimates with 90% Confidence Intervals",
+    x = "Model",
+    y = "Mean Estimate (€/month)",
+    color = "Group"
+  ) +
+  theme_minimal() +
+  theme(
+    legend.position = "right",
+    strip.text.y = element_text(size = 12, face = "bold"),  # Adjust facet labels
+    axis.text.y = element_text(size = 10)  # Ensure readable model names
+  )