diff --git a/Scripts/MAKE_FILE.R b/Scripts/MAKE_FILE.R
index dd201e5a3475b9421e8c00a235ed162501446440..d0f2d3397734b4f515e516ca9d932854fcc2346c 100644
--- a/Scripts/MAKE_FILE.R
+++ b/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")
+
+
+
+
diff --git a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Pred.R b/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Pred.R
index 1f1b6c04a747dfdd0c51238e8d5884e906e78511..b621e7321c7fd6c8bd3fb3c364ea738e28650019 100644
--- a/Scripts/mxl/Prediction models/Split_samples/mxl_wtp_space_caseM_Control_Pred.R
+++ b/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"))
diff --git a/Scripts/ols/ols_consequentiality.R b/Scripts/ols/ols_consequentiality.R
index 71bb7d02450349a21120f741facb7398a86cb15d..981a90149b5a145a26d659591950ae4994dde4b8 100644
--- a/Scripts/ols/ols_consequentiality.R
+++ b/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)
diff --git a/Scripts/propensity_score_matching.R b/Scripts/propensity_score_matching.R
new file mode 100644
index 0000000000000000000000000000000000000000..20c99b88879383e8ed7e368df70892b15f4fd728
--- /dev/null
+++ b/Scripts/propensity_score_matching.R
@@ -0,0 +1,25 @@
+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")
+
diff --git a/Scripts/visualize_distr_spli.R b/Scripts/visualize_distr_spli.R
index a299d8b58309d8fb161da07b134f7f8aecb4f561..ad73703cd1df9ba1dc1b0c3750653105156637a3 100644
--- a/Scripts/visualize_distr_spli.R
+++ b/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
+ )