diff --git a/Scripts/compare_split_samples.R b/Scripts/compare_split_samples.R
index ab93b95f547579cd013578fbe609128faf50e4ab..78c478617ebdef0efc2ac513091e57fa985ac323 100644
--- a/Scripts/compare_split_samples.R
+++ b/Scripts/compare_split_samples.R
@@ -1,203 +1,206 @@
-#### Load WTP models for different treatments #####
-
-# Case A
-
-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")
-
-
-mxl_info_compare <- as.data.frame(mxl_tr$estimate)
-mxl_info_compare[2] <- as.data.frame(mxl_vol_tr$estimate)
-mxl_info_compare[3] <- as.data.frame(mxl_not_tr$estimate)
-
-alpha = 0.1
-mxl_info_compare$margin_of_error1 <- qnorm(1-alpha)*mxl_tr$robse
-mxl_info_compare$margin_of_error2 <- qnorm(1-alpha)*mxl_vol_tr$robse
-mxl_info_compare$margin_of_error3 <- qnorm(1-alpha)*mxl_not_tr$robse
-
-mxl_info_compare <- rownames_to_column(mxl_info_compare, "Coefficent")
-colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
- "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
-
-
-mxl_melt_info <- melt(mxl_info_compare[1:4], id = "Coefficent")
-mxl_melt_info$ME <- mxl_info_compare$Margin_of_error_TR
-mxl_melt_info$ME[9:16] <- mxl_info_compare$Margin_of_error_VOL_TR
-mxl_melt_info$ME[17:24] <- mxl_info_compare$Margin_of_error_NOT_TR
-
-
-# Figure paper compare info treatments #
-ggplot(data=mxl_melt_info, aes(x=Coefficent, y=abs(value), fill=variable)) +
- geom_bar(stat="identity", position='dodge', width = 0.9) +
- geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
- ylab("Absolute Value") +
- xlab("Coefficient") +
- scale_x_discrete(guide = guide_axis(angle = 45)) +
- scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
- theme(legend.position = c(0.85, 0.8))
-
-ggsave("Figures/compare_split_samplesA.png", dpi = "print", width = 7, height = 5)
-
-# Case A log rent
-
-mxl_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Treated A log")
-mxl_vol_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Treated A log")
-mxl_not_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Not_Treated A log")
-
-
-mxl_info_compare <- as.data.frame(mxl_tr$estimate)
-mxl_info_compare[2] <- as.data.frame(mxl_vol_tr$estimate)
-mxl_info_compare[3] <- as.data.frame(mxl_not_tr$estimate)
-
-alpha = 0.05
-mxl_info_compare$margin_of_error1 <- qnorm(1-alpha)*mxl_tr$robse
-mxl_info_compare$margin_of_error2 <- qnorm(1-alpha)*mxl_vol_tr$robse
-mxl_info_compare$margin_of_error3 <- qnorm(1-alpha)*mxl_not_tr$robse
-
-mxl_info_compare <- rownames_to_column(mxl_info_compare, "Coefficent")
-colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
- "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
-
-
-mxl_melt_info_log <- melt(mxl_info_compare[1:4], id = "Coefficent")
-mxl_melt_info_log$ME <- mxl_info_compare$Margin_of_error_TR
-mxl_melt_info_log$ME[8:14] <- mxl_info_compare$Margin_of_error_VOL_TR
-mxl_melt_info_log$ME[15:21] <- mxl_info_compare$Margin_of_error_NOT_TR
-
-mxl_melt_info_log <- mxl_melt_info_log %>% mutate(ME = case_when(abs(value) < 1 ~ ME*100, TRUE~ME),
- value = case_when(abs(value) < 1 ~ value*100, TRUE~value))
-
-# Figure paper compare info treatments #
-ggplot(data=mxl_melt_info_log, aes(x=Coefficent, y=abs(value), fill=variable)) +
- geom_bar(stat="identity", position='dodge', width = 0.9) +
- geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
- ylab("Absolute Value") +
- xlab("Coefficient") +
- scale_x_discrete(guide = guide_axis(angle = 45)) +
- scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
- theme(legend.position = c(0.85, 0.8))
-
-
-
-### Compare clogit models prefspace
-
-clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR")
-clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR")
-clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR")
-
-
-clogit_compare <- as.data.frame(clogit_tr$estimate)
-clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
-clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
-
-alpha = 0.1
-clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
-
-clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
-colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
- "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
-
-
-clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
-clogit_melt$ME <-clogit_compare$Margin_of_error_TR
-clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
-clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
-
-# Scale up rent and wd coefficients for visibility
-clogit_melt <- clogit_melt %>% mutate(ME = case_when(value <= 0.1 ~ ME*10, TRUE ~ ME),
- value = case_when(value <= 0.1 ~ value*10, TRUE ~ value))
-
-# Figure clogit prefspace
-ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
- geom_bar(stat="identity", position='dodge', width = 0.9) +
- geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
- ylab("Absolute Value") +
- xlab("Coefficient") +
- scale_x_discrete(guide = guide_axis(angle = 45)) +
- scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
- theme(legend.position = c(0.85, 0.8))
-
-
-## Clogit WTP space
-
-clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR wtp")
-clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR wtp")
-clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR wtp")
-
-
-clogit_compare <- as.data.frame(clogit_tr$estimate)
-clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
-clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
-
-alpha = 0.1
-clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
-
-clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
-colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
- "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
-
-
-clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
-clogit_melt$ME <-clogit_compare$Margin_of_error_TR
-clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
-clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
-
-# Scale up rent and wd coefficients for visibility
-clogit_melt <- clogit_melt %>% mutate(ME = case_when(abs(value) <= 0.5 ~ ME*1000, abs(value) <= 2 ~ ME*10, TRUE ~ ME),
- value = case_when(abs(value) <= 0.5 ~ abs(value)*1000, value <= 2 ~ value*10, TRUE ~ value))
-
-# Figure clogit prefspace
-ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
- geom_bar(stat="identity", position='dodge', width = 0.9) +
- geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
- ylab("Absolute Value") +
- xlab("Coefficient") +
- scale_x_discrete(guide = guide_axis(angle = 45)) +
- scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
- theme(legend.position = c(0.85, 0.8))
-
-
-## Clogit WTP space log rent
-
-clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR wtp log")
-clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR wtp log")
-clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR wtp log")
-
-
-clogit_compare <- as.data.frame(clogit_tr$estimate)
-clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
-clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
-
-alpha = 0.1
-clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
-
-clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
-colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
- "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
-
-
-clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
-clogit_melt$ME <-clogit_compare$Margin_of_error_TR
-clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
-clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
-
-# Scale up rent and wd coefficients for visibility
-clogit_melt <- clogit_melt %>% mutate(ME = case_when(Coefficent != "mu_rent" ~ ME*1000, TRUE ~ ME),
- value = case_when(Coefficent != "mu_rent" ~ value*1000, TRUE ~ value))
-
-# Figure clogit prefspace
-ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
- geom_bar(stat="identity", position='dodge', width = 0.9) +
- geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
- ylab("Absolute Value") +
- xlab("Coefficient") +
- scale_x_discrete(guide = guide_axis(angle = 45)) +
- scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
- theme(legend.position = c(0.85, 0.8))
+#### Load WTP models for different treatments #####
+
+# Case A
+
+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")
+
+
+mxl_info_compare <- as.data.frame(mxl_tr$estimate)
+mxl_info_compare[2] <- as.data.frame(mxl_vol_tr$estimate)
+mxl_info_compare[3] <- as.data.frame(mxl_not_tr$estimate)
+
+alpha = 0.1
+mxl_info_compare$margin_of_error1 <- qnorm(1-alpha)*mxl_tr$robse
+mxl_info_compare$margin_of_error2 <- qnorm(1-alpha)*mxl_vol_tr$robse
+mxl_info_compare$margin_of_error3 <- qnorm(1-alpha)*mxl_not_tr$robse
+
+mxl_info_compare <- rownames_to_column(mxl_info_compare, "Coefficent")
+colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
+ "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
+
+
+mxl_melt_info <- melt(mxl_info_compare[1:4], id = "Coefficent")
+mxl_melt_info$ME <- mxl_info_compare$Margin_of_error_TR
+
+coeff_num <- nrow(mxl_info_compare)
+
+mxl_melt_info$ME[(coeff_num+1):(coeff_num*2)] <- mxl_info_compare$Margin_of_error_VOL_TR
+mxl_melt_info$ME[(coeff_num*2+1):(coeff_num*3)] <- mxl_info_compare$Margin_of_error_NOT_TR
+
+
+# Figure paper compare info treatments #
+ggplot(data=mxl_melt_info, aes(x=Coefficent, y=abs(value), fill=variable)) +
+ geom_bar(stat="identity", position='dodge', width = 0.9) +
+ geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+ ylab("Absolute Value") +
+ xlab("Coefficient") +
+ scale_x_discrete(guide = guide_axis(angle = 45)) +
+ scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
+ theme(legend.position = c(0.85, 0.8))
+
+ggsave("Figures/compare_split_samplesA.png", dpi = "print", width = 7, height = 5)
+
+# Case A log rent
+
+mxl_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Treated A log")
+mxl_vol_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Vol_Treated A log")
+mxl_not_tr <- apollo_loadModel("Estimation_results/mxl/Split_samples/MXL_wtp Not_Treated A log")
+
+
+mxl_info_compare <- as.data.frame(mxl_tr$estimate)
+mxl_info_compare[2] <- as.data.frame(mxl_vol_tr$estimate)
+mxl_info_compare[3] <- as.data.frame(mxl_not_tr$estimate)
+
+alpha = 0.05
+mxl_info_compare$margin_of_error1 <- qnorm(1-alpha)*mxl_tr$robse
+mxl_info_compare$margin_of_error2 <- qnorm(1-alpha)*mxl_vol_tr$robse
+mxl_info_compare$margin_of_error3 <- qnorm(1-alpha)*mxl_not_tr$robse
+
+mxl_info_compare <- rownames_to_column(mxl_info_compare, "Coefficent")
+colnames(mxl_info_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
+ "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
+
+
+mxl_melt_info_log <- melt(mxl_info_compare[1:4], id = "Coefficent")
+mxl_melt_info_log$ME <- mxl_info_compare$Margin_of_error_TR
+mxl_melt_info_log$ME[8:14] <- mxl_info_compare$Margin_of_error_VOL_TR
+mxl_melt_info_log$ME[15:21] <- mxl_info_compare$Margin_of_error_NOT_TR
+
+mxl_melt_info_log <- mxl_melt_info_log %>% mutate(ME = case_when(abs(value) < 1 ~ ME*100, TRUE~ME),
+ value = case_when(abs(value) < 1 ~ value*100, TRUE~value))
+
+# Figure paper compare info treatments #
+ggplot(data=mxl_melt_info_log, aes(x=Coefficent, y=abs(value), fill=variable)) +
+ geom_bar(stat="identity", position='dodge', width = 0.9) +
+ geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+ ylab("Absolute Value") +
+ xlab("Coefficient") +
+ scale_x_discrete(guide = guide_axis(angle = 45)) +
+ scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
+ theme(legend.position = c(0.85, 0.8))
+
+
+
+### Compare clogit models prefspace
+
+clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR")
+clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR")
+clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR")
+
+
+clogit_compare <- as.data.frame(clogit_tr$estimate)
+clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
+clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
+
+alpha = 0.1
+clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
+clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
+clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
+
+clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
+colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
+ "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
+
+
+clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
+clogit_melt$ME <-clogit_compare$Margin_of_error_TR
+clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
+clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
+
+# Scale up rent and wd coefficients for visibility
+clogit_melt <- clogit_melt %>% mutate(ME = case_when(value <= 0.1 ~ ME*10, TRUE ~ ME),
+ value = case_when(value <= 0.1 ~ value*10, TRUE ~ value))
+
+# Figure clogit prefspace
+ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
+ geom_bar(stat="identity", position='dodge', width = 0.9) +
+ geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+ ylab("Absolute Value") +
+ xlab("Coefficient") +
+ scale_x_discrete(guide = guide_axis(angle = 45)) +
+ scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
+ theme(legend.position = c(0.85, 0.8))
+
+
+## Clogit WTP space
+
+clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR wtp")
+clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR wtp")
+clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR wtp")
+
+
+clogit_compare <- as.data.frame(clogit_tr$estimate)
+clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
+clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
+
+alpha = 0.1
+clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
+clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
+clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
+
+clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
+colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
+ "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
+
+
+clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
+clogit_melt$ME <-clogit_compare$Margin_of_error_TR
+clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
+clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
+
+# Scale up rent and wd coefficients for visibility
+clogit_melt <- clogit_melt %>% mutate(ME = case_when(abs(value) <= 0.5 ~ ME*1000, abs(value) <= 2 ~ ME*10, TRUE ~ ME),
+ value = case_when(abs(value) <= 0.5 ~ abs(value)*1000, value <= 2 ~ value*10, TRUE ~ value))
+
+# Figure clogit prefspace
+ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
+ geom_bar(stat="identity", position='dodge', width = 0.9) +
+ geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+ ylab("Absolute Value") +
+ xlab("Coefficient") +
+ scale_x_discrete(guide = guide_axis(angle = 45)) +
+ scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
+ theme(legend.position = c(0.85, 0.8))
+
+
+## Clogit WTP space log rent
+
+clogit_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A TR wtp log")
+clogit_vol_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A VOL TR wtp log")
+clogit_not_tr <- apollo_loadModel("Estimation_results/clogit/Clogit case A NOT TR wtp log")
+
+
+clogit_compare <- as.data.frame(clogit_tr$estimate)
+clogit_compare[2] <- as.data.frame(clogit_vol_tr$estimate)
+clogit_compare[3] <- as.data.frame(clogit_not_tr$estimate)
+
+alpha = 0.1
+clogit_compare$margin_of_error1 <- qnorm(1-alpha)*clogit_tr$robse
+clogit_compare$margin_of_error2 <- qnorm(1-alpha)*clogit_vol_tr$robse
+clogit_compare$margin_of_error3 <- qnorm(1-alpha)*clogit_not_tr$robse
+
+clogit_compare <- rownames_to_column(clogit_compare, "Coefficent")
+colnames(clogit_compare) <- c("Coefficent", "Estimate_TR", "Estimate_VOL_TR", "Estimate_NOT_TR", "Margin_of_error_TR",
+ "Margin_of_error_VOL_TR", "Margin_of_error_NOT_TR")
+
+
+clogit_melt <- melt(clogit_compare[1:4], id = "Coefficent")
+clogit_melt$ME <-clogit_compare$Margin_of_error_TR
+clogit_melt$ME[5:8] <- clogit_compare$Margin_of_error_VOL_TR
+clogit_melt$ME[9:12] <- clogit_compare$Margin_of_error_NOT_TR
+
+# Scale up rent and wd coefficients for visibility
+clogit_melt <- clogit_melt %>% mutate(ME = case_when(Coefficent != "mu_rent" ~ ME*1000, TRUE ~ ME),
+ value = case_when(Coefficent != "mu_rent" ~ value*1000, TRUE ~ value))
+
+# Figure clogit prefspace
+ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
+ geom_bar(stat="identity", position='dodge', width = 0.9) +
+ geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+ ylab("Absolute Value") +
+ xlab("Coefficient") +
+ scale_x_discrete(guide = guide_axis(angle = 45)) +
+ scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
+ theme(legend.position = c(0.85, 0.8))
diff --git a/Scripts/create_tables.R b/Scripts/create_tables.R
index cceb70ddc0880cf762c5b7e9843e83433514ba4c..a56556e2ca61f157fab01c38fee3242d5b393bd7 100644
--- a/Scripts/create_tables.R
+++ b/Scripts/create_tables.R
@@ -242,6 +242,27 @@ texreg(c(case_B_cols_NR[1], remGOF(case_B_cols_NR[2:6])),
file="Tables/mxl/case_B_rent_INT_NR.tex")
+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")
+
+not_tr <- quicktexregapollo(mxl_not_tr)
+treated <- quicktexregapollo(mxl_tr)
+opt_tr <- quicktexregapollo(mxl_vol_tr)
+
+texreg(c(not_tr, treated, opt_tr),
+ custom.coef.map = list("mu_natural" = "Naturalness", "mu_walking" = "Walking Distance", "mu_rent" = "Rent",
+ "ASC_1" = "ASC 1","ASC_2" = "ASC 2", "_natural" = "Naturalness", "nat" = "Naturalness",
+ "wd" = "Walking Distance", "asc" = "ASC SQ", "ASC" ="ASC SQ",
+ "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance",
+ "sig_natural" = "Naturalness SD", "sig_walking" = "Walking Distance SD", "sig_rent" = "Rent SD",
+ "sig_ASC_1" = "ASC 1 SD", "sig_ASC_2" = "ASC 2 SD"),
+ custom.model.names = c("Not Treated", "Treated", "Optional Treatment"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
+ stars = c(0.01, 0.05, 0.1), float.pos="tb",
+ label = "tab:mxl_C_NR",
+ caption = "Results of mixed logit models for case A.",
+ file="Tables/mxl/case_A_splits.tex")
+
# Main model
# texreg(l=list(mxl_wtp_case_a_rentINT),
# custom.coef.map = list("mu_natural" = "Naturalness", "mu_walking" = "Walking Distance", "mu_rent" = "Rent",
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 e31c769eddf9220d579df0d801a8f3a6c85c0895..c9afe6264f01f6829b1ad706e7fd7f6ab77508e0 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
@@ -49,11 +49,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -63,7 +65,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -77,7 +79,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -99,13 +102,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
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 b621e7321c7fd6c8bd3fb3c364ea738e28650019..825e8c113341672065e33e6ec527ed98ff93ce97 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
@@ -50,11 +50,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -64,7 +66,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -78,7 +80,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -100,13 +103,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
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 033604b4950ba7d6a7dee0f816440731cf6485e5..274feacca1f2c02b9281f719016f231b2ea9f1ed 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
@@ -49,11 +49,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -63,7 +65,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -77,7 +79,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -99,13 +102,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
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 92767fe966abe29e913a6568a83a0899005c0975..31931f6c57d4bad33c0a9e66d05d0c24de4be103 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
@@ -49,11 +49,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -63,7 +65,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -77,7 +79,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -99,13 +102,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
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 9ebe574cff14786ffc6aba841e96fe06a5e04105..9f6b408cbb013a89bbfa495242485942075137fe 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
@@ -49,11 +49,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = -5,
+ ASC_2 = -5,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -63,7 +65,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -77,7 +79,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -99,13 +102,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
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 dc5791d135f41ed5065319a92be2eaa98d76fc45..ac4e59281f94bdbf1da2fe6bf3fc390596b72a65 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
@@ -49,11 +49,13 @@ apollo_control = list(
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -63,7 +65,7 @@ apollo_draws = list(
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -77,7 +79,8 @@ apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -99,13 +102,13 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_NR.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_not_tr_A_NR.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_NR.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_not_tr_A_NR.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_log.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_not_tr_A_log.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_log.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_not_tr_A_log.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_NR.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_tr_A_NR.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_NR.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_tr_A_NR.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_log.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_tr_A_log.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_log.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_tr_A_log.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A_log.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_vol_tr_A_log.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A_log.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_vol_tr_A_log.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat_NR.R b/Scripts/mxl/Split_samples/Old models/mxl_wtp_space_vol_treat_NR.R
similarity index 100%
rename from Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat_NR.R
rename to Scripts/mxl/Split_samples/Old models/mxl_wtp_space_vol_treat_NR.R
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R
index efe460191ae6adf5ae9f007c7c33cf6be6a6caea..934442f262729541de100295ec5ef2f581ae4d24 100644
--- a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A.R
@@ -36,11 +36,13 @@ database <- database_full %>%
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -50,7 +52,7 @@ database <- database_full %>%
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -64,7 +66,8 @@ database <- database_full %>%
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -86,13 +89,13 @@ database <- database_full %>%
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
@@ -129,7 +132,7 @@ database <- database_full %>%
mxl_wtp_not_tr_a = 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/Split_samples/mxl_wtp_space_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R
index d238b0d9348bef854cb5ad771bfa067bac14dd40..2b488e6f1195228adaba40a465a5df0f9cf80946 100644
--- a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A.R
@@ -36,11 +36,13 @@ database <- database_full %>%
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 2)
+ sig_ASC_1 = 2,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -50,7 +52,7 @@ database <- database_full %>%
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -64,7 +66,8 @@ database <- database_full %>%
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -86,13 +89,13 @@ database <- database_full %>%
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
@@ -129,7 +132,7 @@ database <- database_full %>%
mxl_wtp_tr_a = 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/Split_samples/mxl_wtp_space_vol_tr_A.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R
index 52d1642d273b46e01ba416b0bc7f8bba02c093b8..12a67621826c66db7e00065c56f6b0e3ab6cb606 100644
--- a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_tr_A.R
@@ -36,11 +36,13 @@ database <- database_full %>%
apollo_beta=c(mu_natural = 15,
mu_walking = -1,
mu_rent = -2,
- ASC_sq = 0,
+ ASC_1 = 0,
+ ASC_2 = 0,
sig_natural = 15,
sig_walking = 2,
sig_rent = 2,
- sig_ASC_sq = 0)
+ sig_ASC_1 = 0,
+ sig_ASC_2 = 0)
### specify parameters that should be kept fixed, here = none
apollo_fixed = c()
@@ -50,7 +52,7 @@ database <- database_full %>%
interDrawsType = "sobol",
interNDraws = n_draws,
interUnifDraws = c(),
- interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+ interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc1", "draws_asc2"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
@@ -64,7 +66,8 @@ database <- database_full %>%
randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
- randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+ randcoeff[["b_ASC_1"]] = ASC_1 + sig_ASC_1 * draws_asc1
+ randcoeff[["b_ASC_2"]] = ASC_2 + sig_ASC_2 * draws_asc2
return(randcoeff)
}
@@ -86,13 +89,13 @@ database <- database_full %>%
# Define utility functions here:
V = list()
- V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+ V[['alt1']] = -b_mu_rent*(b_ASC_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
Rent_1)
- V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+ V[['alt2']] = -b_mu_rent*(b_ASC_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
Rent_2)
- V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+ V[['alt3']] = -b_mu_rent*(b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
Rent_3)
@@ -129,7 +132,7 @@ database <- database_full %>%
mxl_wtp_vol_tr_a = 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/visualize_distr_spli.R b/Scripts/visualize_distr_spli.R
index ecce6d20436f3637ecb0fbbec1d07ca585869934..9fc1bdf0caa59b119436cd164864ef3804e15e64 100644
--- a/Scripts/visualize_distr_spli.R
+++ b/Scripts/visualize_distr_spli.R
@@ -129,14 +129,13 @@ z_value <- (MXL_wtp_Treated_Pred_model$estimate["mu_natural"] - MXL_wtp_Control_
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))))
+ (sqrt((model_1$robse["mu_natural"]^2 + model_2$robse["mu_natural"]^2)))
return(z_value)
}
# Example of using the function
-z_value <- calculate_z_value(mxl_tr, mxl_vol_tr)
+z_value <- calculate_z_value(MXL_wtp_opt_treated, MXL_wtp_opt_not_treated)
# Print the z-value
z_value
@@ -168,7 +167,7 @@ ggplot(plot_data, aes(x = x, y = density, color = Group, linetype = Prediction))
#### Compare means
# Create an empty data frame
# Function to compute CI dynamically
-calculate_ci <- function(mean, se, confidence = 0.90) {
+calculate_ci <- function(mean, se, confidence = 0.95) {
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
@@ -187,7 +186,7 @@ for (i in seq_along(models)) {
se_robust <- model$robse["mu_natural"]
# Compute 90% confidence intervals dynamically
- ci_values <- calculate_ci(mean_value, se_robust, confidence = 0.90)
+ ci_values <- calculate_ci(mean_value, se_robust, confidence = 0.95)
# Append to data frame
temp_data <- data.frame(
@@ -232,3 +231,185 @@ ggplot(ci_data, aes(x = Model, y = Mean, color = Group)) +
strip.text.y = element_text(size = 12, face = "bold"), # Adjust facet labels
axis.text.y = element_text(size = 10) # Ensure readable model names
)
+
+
+#### Do and plot z tests ####
+
+# Function to calculate Z-test and confidence intervals
+calculate_z_test <- function(model_1, model_2, coef="mu_natural") {
+ mean_diff <- model_1$estimate[coef] - model_2$estimate[coef]
+ se_diff <- sqrt(model_1$robse[coef]^2 + model_2$robse[coef]^2)
+
+ # Compute confidence intervals
+ lower_ci <- mean_diff - 1.96* se_diff
+ upper_ci <- mean_diff + 1.96 * se_diff
+
+ # Compute Z-score
+ z_value <- mean_diff / se_diff
+ p_value <- 2 * (1 - pnorm(abs(z_value))) # Two-tailed p-value
+
+ return(data.frame(
+ Model1 = names(model_1$estimate[coef]),
+ Model2 = names(model_2$estimate[coef]),
+ Mean_Difference = mean_diff,
+ Lower_CI = lower_ci,
+ Upper_CI = upper_ci,
+ Z_Value = z_value,
+ P_Value = p_value
+ ))
+}
+
+# Define groups
+not_predicted_models <- models[names(models) %in% c("Treated Not Pred", "Opt Not Pred", "Control Not Pred")]
+predicted_models <- models[names(models) %in% c("Treated Pred", "Opt Pred", "Control Pred")]
+split_sample_models <- models[names(models) %in% c("Treated Split Sample", "Opt Split Sample", "Control Split Sample")]
+
+
+
+ # Function to perform Z-tests within a group while ensuring the correct order
+ perform_z_tests_within_group <- function(model_group, group_name) {
+ model_names <- names(model_group)
+ z_test_results <- data.frame()
+
+ # Define the correct comparison order
+ comparison_order <- c("Treated", "Optional", "Control") # Order: Treated vs Control, Treated vs Optional, Optional vs Control
+
+ # Reorder the models according to the desired comparison order
+ ordered_models <- model_group[order(match(sub(" .*", "", model_names), comparison_order))]
+
+ # Perform Z-tests in the correct order
+ for (i in seq_along(ordered_models)) {
+ for (j in seq_along(ordered_models)) {
+ if (i < j) { # Avoid duplicate comparisons
+ test_result <- calculate_z_test(ordered_models[[i]], ordered_models[[j]], coef="mu_natural")
+ test_result$Comparison <- paste(names(ordered_models)[i], "vs", names(ordered_models)[j]) # Add model pair label
+ test_result$Group <- group_name # Label the group
+ z_test_results <- bind_rows(z_test_results, test_result)
+ }
+ }
+ }
+ return(z_test_results)
+ }
+
+ # Function to perform cross-group tests while ensuring "Predicted vs Not Predicted"
+ perform_cross_group_tests <- function(not_predicted, predicted) {
+ z_test_results <- data.frame()
+ comparison_order <- c("Treated", "Optional", "Control") # Ensure correct order
+
+ category_map <- list(
+ "Treated Pred" = "Treated Not Pred",
+ "Opt Pred" = "Opt Not Pred",
+ "Control Pred" = "Control Not Pred"
+ )
+
+ # Sort models based on the predefined comparison order
+ ordered_predicted <- predicted[order(match(sub(" .*", "", names(predicted)), comparison_order))]
+ ordered_not_predicted <- not_predicted[order(match(sub(" .*", "", names(not_predicted)), comparison_order))]
+
+ for (cat_pred in names(category_map)) {
+ cat_not_pred <- category_map[[cat_pred]]
+
+ if (cat_pred %in% names(ordered_predicted) & cat_not_pred %in% names(ordered_not_predicted)) {
+ test_result <- calculate_z_test(ordered_predicted[[cat_pred]], ordered_not_predicted[[cat_not_pred]], coef="mu_natural")
+ test_result$Comparison <- paste(cat_pred, "vs", cat_not_pred) # Now "Predicted vs Not Predicted"
+ test_result$Group <- "Predicted vs Not Predicted"
+ z_test_results <- bind_rows(z_test_results, test_result)
+ }
+ }
+
+ return(z_test_results)
+ }
+
+ # Perform Z-tests within each group while enforcing the correct order
+ z_tests_not_predicted <- perform_z_tests_within_group(not_predicted_models, "Not Predicted")
+ z_tests_predicted <- perform_z_tests_within_group(predicted_models, "Predicted")
+ z_tests_split_sample <- perform_z_tests_within_group(split_sample_models, "Split Sample")
+
+ # Perform cross-group tests
+ z_tests_cross_group <- perform_cross_group_tests(not_predicted_models, predicted_models)
+
+ # Combine all Z-test results
+ z_test_results_all <- bind_rows(
+ z_tests_not_predicted,
+ z_tests_predicted,
+ z_tests_split_sample,
+ z_tests_cross_group
+ )
+
+ # Add significance column
+ z_test_results_all$Significant <- ifelse(z_test_results_all$P_Value < 0.05, "Significant", "Not Significant")
+
+ # Ensure correct ordering of groups (Not Predicted → Predicted → Split Sample → Predicted vs Not Predicted)
+ z_test_results_all$Group <- factor(z_test_results_all$Group, levels = c("Not Predicted", "Predicted", "Split Sample", "Predicted vs Not Predicted"))
+
+ # Separate groups for visualization
+ z_test_results_within <- z_test_results_all %>% filter(Group %in% c("Not Predicted", "Predicted"))
+ z_test_results_splits <- z_test_results_all %>% filter(Group %in% c("Split Sample"))
+ z_test_results_cross <- z_test_results_all %>% filter(Group == "Predicted vs Not Predicted")
+
+ ### **Visualization 1: Within-Group Comparisons (Ordered by Group)** ###
+ ggplot(z_test_results_within, aes(x = reorder(Comparison, Group), y = Mean_Difference, color = Significant)) +
+ geom_point(size = 3) +
+ geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2) +
+ scale_color_manual(values = c("Significant" = "red", "Not Significant" = "blue")) +
+ geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Line at zero
+ coord_flip() + # Flip for readability
+ labs(
+ #title = "Within-Group Z-Test Results for Mean Differences",
+ x = "Model Comparisons",
+ y = "Mean Difference (Naturalness)",
+ color = "Significance"
+ ) +
+ facet_wrap(~ Group, scales = "free_y", ncol = 1) + # Separate groups visually
+ theme_minimal() +
+ theme(legend.position = "right")
+
+ ggsave("Figures/z_test_predictions_between.png", width=8, height=5, dpi="print")
+
+ #### Only split samples
+
+ z_test_results_splits <- z_test_results_splits %>%
+ mutate(Comparison = case_when(
+ grepl("Treated Split Sample vs Control Split Sample", Comparison) ~ "Treated vs Control",
+ grepl("Treated Split Sample vs Opt Split Sample", Comparison) ~ "Treated vs Optional",
+ grepl("Control Split Sample vs Opt Split Sample", Comparison) ~ "Control vs Optional",
+ TRUE ~ Comparison # Keep original if not matched
+ ))
+
+ ggplot(z_test_results_splits, aes(x = reorder(Comparison, Group), y = Mean_Difference, color = Significant)) +
+ geom_point(size = 3) +
+ geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2) +
+ scale_color_manual(values = c("Significant" = "red", "Not Significant" = "blue")) +
+ geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Line at zero
+ coord_flip() + # Flip for readability
+ labs(
+ # title = "Within-Group Z-Test Results for Mean Differences",
+ x = "Model Comparisons",
+ y = "Mean Difference (Naturalness)",
+ color = "Significance"
+ ) +
+ theme_minimal() +
+ theme(legend.position = "right")
+
+ ggsave("Figures/z_test_split_samples.png", width=8, height=5, dpi="print")
+
+
+ ### **Visualization 2: Predicted vs Not Predicted Comparisons** ###
+ ggplot(z_test_results_cross, aes(x = Comparison, y = Mean_Difference, color = Significant)) +
+ geom_point(size = 3) +
+ geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = 0.2) +
+ scale_color_manual(values = c("Significant" = "red", "Not Significant" = "blue")) +
+ geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Line at zero
+ coord_flip() + # Flip for readability
+ labs(
+ #title = "Predicted vs Not Predicted (Same Group) Z-Test Results",
+ x = "Model Comparisons",
+ y = "Mean Difference (Naturalness)",
+ color = "Significance"
+ ) +
+ theme_minimal() +
+ theme(legend.position = "right")
+
+ ggsave("Figures/z_test_predictions.png", width=8, height=5, dpi="print")
+
+
\ No newline at end of file