diff --git a/Scripts/compare_split_samples.R b/Scripts/compare_split_samples.R
index 011aa550135ddf2750e2c30e59fca42acda8d3b4..9fbb6064d2547d33c12581360ba840e111a2fbaf 100644
--- a/Scripts/compare_split_samples.R
+++ b/Scripts/compare_split_samples.R
@@ -12,9 +12,9 @@ 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/2)*mxl_tr$robse
-mxl_info_compare$margin_of_error2 <- qnorm(1-alpha/2)*mxl_vol_tr$robse
-mxl_info_compare$margin_of_error3 <- qnorm(1-alpha/2)*mxl_not_tr$robse
+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",
@@ -34,7 +34,7 @@ ggplot(data=mxl_melt_info, aes(x=Coefficent, y=abs(value), fill=variable)) +
   ylab("Absolute Value") +
   xlab("Coefficient") +
   scale_x_discrete(guide = guide_axis(angle = 45)) +
-  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  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)
@@ -51,9 +51,9 @@ 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/2)*mxl_tr$robse
-mxl_info_compare$margin_of_error2 <- qnorm(1-alpha/2)*mxl_vol_tr$robse
-mxl_info_compare$margin_of_error3 <- qnorm(1-alpha/2)*mxl_not_tr$robse
+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",
@@ -75,7 +75,7 @@ ggplot(data=mxl_melt_info_log, aes(x=Coefficent, y=abs(value), fill=variable)) +
   ylab("Absolute Value") +
   xlab("Coefficient") +
   scale_x_discrete(guide = guide_axis(angle = 45)) +
-  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
   theme(legend.position = c(0.85, 0.8)) 
 
 
@@ -92,9 +92,9 @@ 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/2)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha/2)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha/2)*clogit_not_tr$robse
+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",
@@ -117,7 +117,7 @@ ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
   ylab("Absolute Value") +
   xlab("Coefficient") +
   scale_x_discrete(guide = guide_axis(angle = 45)) +
-  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
   theme(legend.position = c(0.85, 0.8)) 
 
 
@@ -133,9 +133,9 @@ 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/2)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha/2)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha/2)*clogit_not_tr$robse
+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",
@@ -158,7 +158,7 @@ ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
   ylab("Absolute Value") +
   xlab("Coefficient") +
   scale_x_discrete(guide = guide_axis(angle = 45)) +
-  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Voluntary Treatment", "Not Treated"), name="Treatment") +
   theme(legend.position = c(0.85, 0.8)) 
 
 
@@ -174,9 +174,9 @@ 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/2)*clogit_tr$robse
-clogit_compare$margin_of_error2 <- qnorm(1-alpha/2)*clogit_vol_tr$robse
-clogit_compare$margin_of_error3 <- qnorm(1-alpha/2)*clogit_not_tr$robse
+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",
@@ -199,5 +199,5 @@ ggplot(data=clogit_melt, aes(x=Coefficent, y=abs(value), fill=variable)) +
   ylab("Absolute Value") +
   xlab("Coefficient") +
   scale_x_discrete(guide = guide_axis(angle = 45)) +
-  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  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/interaction_plots_presi.R b/Scripts/interaction_plots_presi.R
new file mode 100644
index 0000000000000000000000000000000000000000..988847dde811ddf4c952e774cddf8237d73cf429
--- /dev/null
+++ b/Scripts/interaction_plots_presi.R
@@ -0,0 +1,82 @@
+# Create Interaction Term Plot for Presentation
+
+
+alpha = 0.1
+z_value <- qnorm(1-alpha/2)
+
+
+### Interview time A
+plot_interview <- summary(ols_time_spent_control_A)
+plot_interview <- as.data.frame(plot_interview$coefficients)
+plot_interview$ME <- z_value*plot_interview$`Std. Error`
+plot_interview <- rownames_to_column(plot_interview, "Coefficient")
+
+plot_interview <- plot_interview %>% filter(str_detect(Coefficient, "Treatment"))
+
+plot_interview$Coefficient <- c("Treated", "Voluntary Treated")
+
+plot_interview_A <- ggplot(data=plot_interview) +
+  geom_bar(aes(x=Coefficient, y=Estimate, fill=Coefficient), stat="identity",  position='dodge', width = 0.5, alpha=0.7) +
+  geom_errorbar(aes(x=Coefficient, ymin=Estimate-ME, ymax=Estimate+ME), width=0.3, position=position_dodge(0.8)) +
+  scale_x_discrete(guide = guide_axis(angle = 45)) +
+  guides(fill = "none") +
+  xlab("Treatment Group")
+
+### Interview time C
+
+plot_interview <- summary(ols_time_spent_control_C)
+plot_interview <- as.data.frame(plot_interview$coefficients)
+plot_interview$ME <- z_value*plot_interview$`Std. Error`
+plot_interview <- rownames_to_column(plot_interview, "Coefficient")
+
+plot_interview <- plot_interview %>% filter(str_detect(Coefficient, "Treatment"))
+
+plot_interview$Coefficient <- c("No Info 2", "Text 1", "Text 2", "Video 1", "Video 2")
+
+plot_interview_C <- ggplot(data=plot_interview) +
+  geom_bar(aes(x=Coefficient, y=Estimate, fill=Coefficient), stat="identity",  position='dodge', width = 0.5, alpha=0.7) +
+  geom_errorbar(aes(x=Coefficient, ymin=Estimate-ME, ymax=Estimate+ME), width=0.3, position=position_dodge(0.8)) +
+  scale_x_discrete(guide = guide_axis(angle = 45)) +
+  guides(fill = "none") +
+  xlab("Treatment Group")
+
+
+####
+create_interaction_term_plot <- function(ols_summary, treatment_labels, unit) {
+  alpha <- 0.1
+  z_value <- qnorm(1 - alpha / 2)
+  
+  plot_data <- summary(ols_summary)
+  plot_data <- as.data.frame(plot_data$coefficients)
+  plot_data$ME <- z_value * plot_data$`Std. Error`
+  plot_data <- rownames_to_column(plot_data, "Coefficient")
+  
+  plot_data <- plot_data %>% filter(str_detect(Coefficient, "Treatment"))
+  
+  plot_data$Coefficient <- treatment_labels
+  
+  plot <- ggplot(data = plot_data) +
+    geom_bar(aes(x = Coefficient, y = Estimate, fill = Coefficient), stat = "identity", position = 'dodge', width = 0.5, alpha = 0.7) +
+    geom_errorbar(aes(x = Coefficient, ymin = Estimate - ME, ymax = Estimate + ME), width = 0.3, position = position_dodge(0.8)) +
+    scale_x_discrete(guide = guide_axis(angle = 45)) +
+    guides(fill = "none") +
+    xlab("Treatment Group") +
+    ylab(paste0(unit))
+  
+  return(plot)
+}
+
+
+case_A_labels <- c("Treated", "Voluntary Treated")
+case_C_labels <- c("No Info 2", "Text 1", "Text 2", "Video 1", "Video 2")
+
+plot_interview_A <- create_interaction_term_plot(ols_time_spent_control_A, case_A_labels, "Interview Time in Seconds")
+plot_interview_C <- create_interaction_term_plot(ols_time_spent_control_C, case_C_labels, "Interview Time in Seconds")
+
+plot_cc_A <- create_interaction_term_plot(ols_time_cc_control_A, case_A_labels, "Mean Choice Card Time in Seconds")
+plot_cc_C <- create_interaction_term_plot(ols_time_cc_control_C, case_C_labels, "Mean Choice Card Time in Seconds")
+
+plot_mani_A <- create_interaction_term_plot(ols_percentage_correct_control_A, case_A_labels, 
+                                            "Percentage of Correct Quiz Statements")
+plot_mani_C <- create_interaction_term_plot(ols_percentage_correct_control_C, case_C_labels, 
+                                            "Percentage of Correct Quiz Statements")
diff --git a/Scripts/nat_interaction_plot.R b/Scripts/nat_interaction_plot.R
new file mode 100644
index 0000000000000000000000000000000000000000..b29594ac72624b73fd5804d340fa5897842ad858
--- /dev/null
+++ b/Scripts/nat_interaction_plot.R
@@ -0,0 +1,26 @@
+df_NR_C <- as.data.frame(mxl_wtp_NR_case_c_rentINT$estimate)
+alpha = 0.1
+z_value <- qnorm(1-alpha)
+df_NR_C$margin_of_error1 <- z_value*mxl_wtp_NR_case_c_rentINT$robse
+
+df_NR_C <- rownames_to_column(df_NR_C, "Coefficent")
+
+colnames(df_NR_C) <- c("Coefficent", "Estimate", "Margin_of_error")
+
+df_NR_C_melt <- melt(df_NR_C[1:2], id = "Coefficent")
+
+df_NR_C_melt$ME <- df_NR_C$Margin_of_error
+
+df_NR_C_melt <- df_NR_C_melt %>%
+  filter(str_starts(Coefficent, "mu_nat_"))
+
+df_NR_C_melt$Coefficent <- c("NR", "Video 1", "Video 2", "No Info 2", "Text 1", "Text 2")
+
+ggplot(data=df_NR_C_melt, 
+       aes(x=factor(Coefficent, levels = c("NR", "No Info 2", "Text 1", "Text 2", "Video 1", "Video 2")),
+          y=value, fill=variable)) +
+  geom_bar(aes(fill=Coefficent), stat="identity",  position='dodge', width = 0.9, alpha=0.7) +
+  geom_errorbar(aes(x=Coefficent, ymin=value-ME, ymax=value+ME), width=0.3, position=position_dodge(0.8)) +
+  ylab("Interaction Effect WTP Naturalness") +
+  xlab("Interaction") +
+  scale_x_discrete(guide = guide_axis(angle = 45))