Tax_compliance_analysis.Rmd

---
title: "Increasing Property Tax Payment Rates on Time Through Simplification of Communication with Citizens"
author: "Sloboda - Pavlovský - Sičáková-Beblavá - Szabo"
date: "2024-07-05"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Read Me

This file replicates the analysis of the following paper: "*Increasing Property Tax Payment Rates on Time Through Simplification of Communication with Citizens*" submitted to [Journal of Economic Behavior & Organization](https://www.sciencedirect.com/journal/journal-of-economic-behavior-and-organization). The file replicates the results shown in the main text as well as in the appendix.

## Dataset Preparation

```{r message=FALSE}
rm(list = ls())

setwd("/Users/jakubszabo/Library/CloudStorage/OneDrive-Personal/Desktop/FSEV/FSEV_PC/Tax_compliance")

# Required Libraries 
library(tidyverse)
library(rio)
library(estimatr) # lm_robust 
library(margins)
library(modelsummary)
library(sjPlot)
library(gridExtra)

# Importing Dataset
raw_data <- import("tax_compliance.csv")

# Renaming variable names from Slovak to English 
data <- raw_data %>% 
  dplyr::rename(
    ID = REFERENCNECISLO, 
    Gender = POHLAVIE, 
    Age = VEK,
    Tax_total = DANCELKOM,
    Tax_paid = UHRADENACIASTKA, 
    Tax_paid_before_duedate = UHRADENACIASTKADOSPLATNOSTI,
    Tax_paid_after_duedate = UHRADENACIASTKAPOSPLATNOSTI, 
    Tax_outstanding = ZOSTAVAUHRADIT, 
    Treatment = INTERVENCIA,
    Tax_declaration_21 = PRIZNANIE2021A,
    Admin_ID = CISLOSPRAVKYNE
    ) 

# Recoding Gender Variable 
data <- data %>% 
  mutate(Gender = case_match(Gender, 
    "zena" ~ "F", 
    "muz" ~ "M"
  )) %>% 
  mutate(sex = if_else(Gender == "M", 0, 1)) %>% 
  relocate(sex, .after = Gender)

data$Gender <- factor(data$Gender, levels = c("M", "F"))

# Recoding Treatment Variable 
data <- data %>%  
 mutate(Treatment = case_match(Treatment,
                                "Export DZN A" ~ "A",
                                "Export DZN B" ~ "B",
                                "Export DZN C" ~ "C",
                                "Export DZN D" ~ "D",
                                "Export DZN E" ~ "E")) %>% 
  mutate(Treatment_full = case_match(Treatment,
                                "A" ~ "Control",
                                "B" ~ "Simplification",
                                "C" ~ "Simplification + Social Norm",
                                "D" ~ "Simplification + Deterrence",
                                "E" ~ "Simplification + Public Good Appeal")) %>% 
  relocate(Treatment_full, .after = Treatment)


data$Treatment_full <- factor(data$Treatment_full, levels = c("Control", "Simplification", "Simplification + Social Norm", "Simplification + Deterrence", "Simplification + Public Good Appeal"))

# Setting as Factor 
data <- data %>% 
  mutate(Admin_ID = as.factor(Admin_ID)) %>% 
  mutate(Tax_declaration_21 = as.factor(Tax_declaration_21)) 

# Creating the Outcome Variable(s) Tax_paid_singlepay & Tax_paid_install 
data <- data %>% 
  mutate(Residul_tax_before_duedate = Tax_total - Tax_paid_before_duedate) %>% 
  mutate(Tax_paid_singlepay = if_else(Residul_tax_before_duedate == 0, 1, 0)) %>% 
   mutate(Tax_paid_late = if_else(Tax_paid_after_duedate == 0, 0, 1)) %>% 
   mutate(Tax_paid_install = if_else(Tax_paid_before_duedate == 0, 0, 1))
```

## Administrator Distribution and Treatment Randomization Check (Appendix B.1)

```{r}
# Materials and procedures - Number of Obs. per Each Treatment for Table 1 
obs_per_treatment <- data %>% 
  group_by(Treatment_full) %>% 
  summarise(n = n())

# print(obs_per_treatment)
```


```{r}
## Table B.1a: Number of Taxpayers per Administrator and Treatment
admin_check <- data %>% 
  group_by(Admin_ID, Treatment) %>% 
  dplyr::summarise(n = n()) %>% 
  pivot_wider(names_from = Treatment, values_from = n) %>% 
  mutate(across(where(is.integer), ~replace_na(., 0))) %>% 
  relocate(A, .before = B) %>% 
  relocate(D, .before = E)

print(admin_check)
```

```{r}
# Filtering wrongly assigned treatments 
admin_check$Treatment <- apply(admin_check, 1, function(row) {
  colnames(admin_check)[which.max(row)]
})

data_subset <- data %>%
  inner_join(admin_check, by = c("Admin_ID", "Treatment")) %>% 
  select(-c(A, B, C, D, E))

# nrow(data) - nrow(data_subset)  

# n_admin <- data_subset %>% 
#  group_by(Admin_ID) %>% 
#  summarise(n = n()) 

# range(n_admin$n)

# data_subset %>% 
#  group_by(Treatment) %>% 
#  summarise(Count = n_distinct(Admin_ID))
```

```{r}
# Table B.1b: Descriptive Statistics of the Selected Variable per Each Treatment
treatment_randomization_check <- data_subset %>% 
  group_by(Treatment) %>% 
  dplyr::summarise(
    `Mean tax total` = round(mean(Tax_total), 2), 
    `Mean tax paid` = round(mean(Tax_paid), 2), 
    `Mean tax outstanding` = round(mean(Tax_outstanding), 2), 
    `Share of people filling tax declaration in 2021 (%)` = round(mean(as.integer(Tax_declaration_21))*100 - 100, 2), 
    `Share of female (%)` = round(mean(sex)*100, 2),
    `Mean age` = round(mean(Age), 2)) %>% 
   pivot_longer(
    cols = -Treatment,
    names_to = "Statistics",
    values_to = "Value"
  ) %>% 
  pivot_wider(
    names_from = Treatment,
    values_from = Value
  )

print(treatment_randomization_check)
```

```{r}
# Table B.1c: Treatment Randomization Check - Total Tax Burden

# anova1 <- aov(Tax_total ~ Treatment, data = data_subset)
# summary(anova1)
# anova1.residuals <- residuals( object = anova1 )   
# qqnorm(anova1.residuals)
# nortest::ad.test(anova1.residuals)

anova_tax <- kruskal.test(Tax_total ~ Treatment, data = data_subset)

anova_tax_df <- data.frame(
  Statistic = anova_tax$statistic,
  `p value` = anova_tax$p.value,
  df = anova_tax$parameter
)

pairwise_wilcox <- pairwise.wilcox.test(
  data_subset$Tax_total, 
  data_subset$Treatment,
  p.adjust.method = "bonferroni"
)

pairwise_wilcox_df <- as.data.frame(pairwise_wilcox$p.value)%>% 
   mutate(across(where(is.numeric), ~ round(., 3)))

pairwise_wilcox_df <- pairwise_wilcox_df %>%
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), "-", format(round(., 2), nsmall = 2)))) %>%
  mutate(across(everything(), ~ ifelse(. == "1.00", "1.00", .)))

print(anova_tax_df)
print(pairwise_wilcox_df)

## Table B.1d: Treatment Randomization Check - Gender 
# chisq_gender <- CrossTable(data_subset$Gender, data_subset$Treatment, 
#            prop.r = FALSE,
#            prop.c = FALSE,
#            prop.t = FALSE,
#            prop.chisq = FALSE,
#            chisq = TRUE)

Gender_Treatment <- table(data_subset$Gender, data_subset$Treatment)
chisq_df <- chisq.test(Gender_Treatment)

chisq_gender_df <- data.frame(
  Statistic = chisq_df$statistic,
  `p value` = chisq_df$p.value,
  df = chisq_df$parameter
)

print(chisq_gender_df)

## Table B.1e: Treatment Randomization Check - Age
# anova2 <- aov(Age ~ Treatment, data = data_subset)
# summary(anova2)
# anova2.residuals <- residuals(object = anova2)   
# qqnorm(anova2.residuals)
# nortest::ad.test(anova2.residuals)

anova_age <- kruskal.test(Age ~ Treatment, data = data_subset)

anova_age_df <- data.frame(
  Statistic = anova_age$statistic,
  `p value` = anova_age$p.value,
  df = anova_age$parameter
)


pairwise_wilcox_age <- pairwise.wilcox.test(
  data_subset$Age, 
  data_subset$Treatment,
  p.adjust.method = "bonferroni"
)

pairwise_wilcox_age_df <- as.data.frame(pairwise_wilcox_age$p.value)%>% 
   mutate(across(where(is.numeric), ~ round(., 3)))

pairwise_wilcox_age_df <- pairwise_wilcox_age_df %>%
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), "-", format(round(., 2), nsmall = 2)))) %>%
  mutate(across(everything(), ~ ifelse(. == "1.00", "1.00", .)))

print(anova_age_df)
print(pairwise_wilcox_age_df)
```

## B.2 Outlier Analysis (Appendix B.2)
```{r}
## Figure B.2a: Share of Taxpayers by Age Decil 
age_mean <- mean(data_subset$Age)

plot_age <- ggplot(data = data_subset, aes(x = Age)) + 
  geom_histogram(fill = "grey", colour = "black", bins = 10)+
  stat_function(fun = function(x) dnorm(x, mean = mean(data$Age), sd = sd(data$Age)),
                color = "black", linetype = "dotted") + 
  geom_vline(xintercept = age_mean, linetype = "dashed", color = "black")+
  theme(
    text = element_text(face = "plain")  
  ) +
  theme_classic() +
  labs(y = "No. of observations", x = "Age") 

print(plot_age)
```

```{r}
# Filtering out age outliers 
data_subset <- data_subset %>% 
  mutate(Age_outlier = if_else(Age < 18 | Age > 100, 1, 0)) 

# data_subset %>% 
#  group_by(Age_outlier) %>% 
#  summarise(n = n())

data_subset <- data_subset %>%
  filter(Age_outlier != 1)
```

```{r}
## Figure B.2b: Distribution of the Tax-Related Variables  

plot_Tax_total <- ggplot(data = data_subset, aes(x = Tax_total)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Total tax")

plot_Tax_total_paid <- ggplot(data = data_subset, aes(x = Tax_paid)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Total tax paid")

plot_Tax_paid_before_duedate  <- ggplot(data = data_subset, aes(x = Tax_paid_before_duedate)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Tax paid before the due date")
  
plot_Tax_paid_after_duedate  <- ggplot(data = data_subset, aes(x = Tax_paid_after_duedate)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Tax paid after the due date")

plot_Tax_outstanding  <- ggplot(data = data_subset, aes(x = Tax_outstanding)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Outstanding tax")

grid.arrange(plot_Tax_total,plot_Tax_total_paid,  plot_Tax_paid_before_duedate, plot_Tax_paid_after_duedate, plot_Tax_outstanding)
```

```{r}
## Filtering out Tax Outliers (+- 2 sd)
sd <- sd(data_subset$Tax_total)
lower_bound2 <- mean(data_subset$Tax_total) - (2 * sd)
upper_bound2 <- mean(data_subset$Tax_total) + (2 * sd)

data_subset <- data_subset %>% 
  mutate(sd2 = if_else(Tax_total > upper_bound2 | Tax_total < lower_bound2, 1, 0))

data_subset_filtered <- data_subset %>% 
  filter(sd2 != 1)
```

```{r}
## Figure B.2c: Distribution of the Tax-Related Variables Without the Outliers 

plot_Tax_total2 <- ggplot(data = data_subset_filtered, aes(x = Tax_total)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Total tax")

plot_Tax_total_paid2 <- ggplot(data = data_subset_filtered, aes(x = Tax_paid)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Total tax paid")

plot_Tax_paid_before_duedate2  <- ggplot(data = data_subset_filtered, aes(x = Tax_paid_before_duedate)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Tax paid before the due date")
  
plot_Tax_paid_after_duedate2  <- ggplot(data = data_subset_filtered, aes(x = Tax_paid_after_duedate)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Tax paid after the due date")

plot_Tax_outstanding2  <- ggplot(data = data_subset_filtered, aes(x = Tax_outstanding)) + 
  geom_density(fill = "grey")+
  theme_classic()+
  labs(y = "",
       x = "Outstanding tax")

grid.arrange(plot_Tax_total2,plot_Tax_total_paid2,  plot_Tax_paid_before_duedate2, plot_Tax_paid_after_duedate2, plot_Tax_outstanding2)
```

## Descriptives Statistics per Treatment after the Cleaning Process (Table 2)
```{r}
table2 <- data_subset_filtered %>% 
  group_by(Treatment_full) %>% 
  summarise(N = n(),
            `Sex (% Female)` = round(mean(sex)*100,2),
            `Mean Age` = round(mean(Age),2),
            `SD Age` = round(sd(Age), 2),
            `Mean Tax` = round(mean(Tax_total),2),
            `SD Tax` = round(sd(Tax_total), 2),
            `Filling Tax Declaration in 2021 (%)` = round(mean(as.numeric(Tax_declaration_21))*100 - 100, 2))

print(table2)
```

## Descriptives Statistics after the Cleaning Process (Appendix B.3)
```{r}
## Table B.3a: Descriptive Statistics of the Numeric Variables (Ungrouped by Treatment) 
desc_stat <- data_subset_filtered %>% 
  select(where(is.numeric)) %>% 
  select(-c(ID, Age_outlier, sd2)) %>%
  psych::describe() %>% 
  select(-vars) 

rownames(desc_stat) <- c("Gender (Female = 1)", 
                         "Age",
                         "Tax total", 
                         "Tax paid", 
                         "Tax paid before duedate",
                         "Tax paid after duedate",
                         "Tax outstanding",
                         "Tax residual before duedate",
                         "Tax paid in a single payment",
                         "Tax paid late",
                         "Tax paid in a single payment / installments"
                         )

desc_stat <- desc_stat %>% 
  select(-c(n, trimmed, mad)) %>%
           mutate(across(where(is.numeric), ~ format(round(., 2), nsmall = 2))) 

print(desc_stat)
```


```{r}
## Figure B.3a: Distribution of Total Tax Burden Across the Control and Treatment Groups
tax_dist_treat <- data_subset_filtered %>% 
  ggplot(aes(x = Treatment_full, y = Tax_total, fill = Treatment_full)) +
  geom_boxplot()+
  geom_violin(alpha = .2)+ 
  stat_summary(geom = "point", fun = "mean", shape = "cross", size = 2)+ 
  # geom_jitter(alpha = .2)+
  theme_classic()+
  labs(y = "Tax total (EUR)", x = "Treatment")+
  theme(legend.position = "none")+
  scale_fill_grey()+
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8))

print(tax_dist_treat)
```

## Analysis
### Average Treatment Effects (DV = Tax paid in a Single Payment or Installments)
```{r}
## Figure 1. Treatment Effects Relative to the Control Group  
# Basic Model with robust SE
model1 <- lm_robust(Tax_paid_install ~ Treatment_full, 
                    data = data_subset_filtered, se_type = "HC3")

# Model with covariates and robust SE
model2 <- lm_robust(Tax_paid_install ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_summary1 <- tidy(model1, conf.int = TRUE)
model_summary1$model <- "Model1"

model_summary2 <- tidy(model2, conf.int = TRUE)
model_summary2$model <- "Model2"

model_all <- rbind(model_summary1, model_summary2) 

model_all <- model_all %>%
  mutate(
    Treatment = gsub("Treatment_full", "", term)
  ) %>% 
  relocate(Treatment, .after = term) 

model_all <- model_all %>% 
  filter(!grepl("Admin|Age|Gender|Tax", term)) %>% 
  mutate(Treatment = if_else(Treatment == "(Intercept)", "Control", Treatment)) %>% 
  select(-term)

model_all$Treatment <- factor(model_all$Treatment, levels = 
                                c("Simplification + Public Good Appeal",
                                   "Simplification + Deterrence", 
                                   "Simplification + Social Norm", 
                                  "Simplification", 
                                  "Control"
                                  ))

figure1 <- model_all %>% 
  filter(!Treatment == "Control") %>% 
  mutate(estimate = estimate*100,
         conf.low = conf.low*100,
         conf.high = conf.high*100) %>% 
  ggplot(aes(y = Treatment, x = estimate, colour = model)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width=0, position= position_dodge(width=0.5)) +
   geom_point(aes(x = estimate, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 2.5), breaks = round(seq(-1, 2.5, by = .5), 2))+
  geom_vline(xintercept = 0, linetype = "dashed")+
  labs(x = expression(atop("Predicted Treatment Effect (95% CI)", italic("Property Tax Payment (single payment or instalments)"), 
      )),  y = "") +
  theme_classic() +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates"), start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(figure1)
# ggsave(figure1, width = 20, height = 12, units = c("cm"), file ="figure1.jpg")
```

```{r}
## Figure C.1a: Predicted Probability of Tax Compliance across the Control and Treatment Groups (95% CI)
pred_model_1 <- prediction(model1, data_subset_filtered)
pred_model_2 <- prediction(model2, data_subset_filtered)

sum_model1 <- plyr::ddply(pred_model_1, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model1$model <- "Model1"

sum_model2 <- plyr::ddply(pred_model_2, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model2$model <- "Model2"

sum_all <- rbind(sum_model1, sum_model2)

plot_pred_prob <- ggplot(sum_all, aes(Treatment_full, mean, colour = model)) +
  geom_errorbar(aes(ymin=mean - (1.96 * mean_se), 
                    ymax= mean + (1.96 * mean_se)), width=0, position= position_dodge(width=0.5)) + 
  geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_y_continuous(limits=c(.92, .95))+ 
  ylab("Predicted Probability of Tax Compliance (95% CI)") +
  xlab("")+ 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", 
                                                  "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", 
                                                  "Covariates"), 
                    start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(plot_pred_prob)
```

### Heterogeneous Treatment Effects (DV = Tax paid on-time in full or in part)
#### Heterogeneous Treatment Effects - Gender
```{r}
## Figure 2. Heterogeneous Treatment Effects by Gender 
model_HTE_gender <- lm_robust(Tax_paid_install ~ Treatment_full * Gender + Tax_total + Tax_declaration_21 + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_gender_inc <- margins(model_HTE_gender, 
                         variables = "Treatment_full", 
                         at = list(Gender = c("M", "F")))

model_HTE_gender_inc <- as.data.frame(summary(model_HTE_gender_inc))

model_HTE_gender_inc$Treatment <- model_HTE_gender_inc$factor
model_HTE_gender_inc$upper95 <- model_HTE_gender_inc$AME + (model_HTE_gender_inc$SE*1.96)
model_HTE_gender_inc$lower95 <- model_HTE_gender_inc$AME - (model_HTE_gender_inc$SE*1.96)

model_HTE_gender_inc <- model_HTE_gender_inc %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    )) %>% 
  mutate(Gender = case_when(
    Gender == "F" ~ "Female", 
    Gender == "M" ~ "Male",
    TRUE ~ Gender))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_gender <- model_HTE_gender_inc %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100) %>% 
  ggplot(aes(colour = Gender)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Gender), position = position_dodge(width=0.5), size=2) +
  theme_classic() + 
  scale_color_grey("Gender") + 
  scale_fill_grey("Gender") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Gender") +
  labs(x = expression(atop("Predicted Treatment Effect (95% CI)", italic("Property Tax Payment (single payment or instalments)"), 
      )),  y = "") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  scale_x_continuous(limits=c(-1, 3), breaks = round(seq(-1, 3, by = .5), 2))

print(plot_inter_gender)
# ggsave(plot_inter_gender, width = 20, height = 12, units = c("cm"), file ="plot_inter_gender.jpg")
```

```{r}
## Figure C.1b: Predicted Probability of Tax Compliance across Treatment Groups and Gender
plot_int_gender <- plot_model(model_HTE_gender, 
           type = "pred",
           terms = c("Treatment_full[all]", "Gender[all]"),
           #colors = "bw",
           legend.title = "Gender")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("Male", "Female"))

print(plot_int_gender)
```

#### Heterogeneous Treatment Effects - Age
```{r}
## Figure 3. Heterogeneous Treatment Effects by Age  
model_HTE_age <- lm_robust(Tax_paid_install ~ Treatment_full * Age + Treatment_full * I(Age^2) + Tax_total + Tax_declaration_21 + Gender, 
                           data = data_subset_filtered, 
                           se_type = "HC3")

model_HTE_age_inc <- margins(model_HTE_age, 
                         variables = "Treatment_full", 
                         at = list(Age = seq(min(data_subset_filtered$Age, na.rm = TRUE), 
                                             max(data_subset_filtered$Age, na.rm = TRUE), 
                                             by = 1)))

model_HTE_age_inc <- as.data.frame(summary(model_HTE_age_inc))

model_HTE_age_inc$upper95 <- model_HTE_age_inc$AME + (model_HTE_age_inc$SE * 1.96)
model_HTE_age_inc$lower95 <- model_HTE_age_inc$AME - (model_HTE_age_inc$SE * 1.96)

model_HTE_age_inc <- model_HTE_age_inc %>%  
  mutate(Treatment = case_when(
    factor == "Treatment_fullSimplification" ~ "Simplification", 
    factor == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    factor == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    factor == "Treatment_fullSimplification + Social Norm" ~ "Simplification + Social Norm",
    TRUE ~ factor 
  ))

model_HTE_age_inc$Treatment <- factor(model_HTE_age_inc$Treatment, levels = c(
                     "Simplification",
                     "Simplification + Social Norm",
                     "Simplification + Deterrence",
                      "Simplification + Public Good Appeal"
                     ))

plot_inter_age_cont <- model_HTE_age_inc %>% 
  mutate(AME = AME * 100,
         lower95 = lower95 * 100,
         upper95 = upper95 * 100) %>%
  ggplot(aes(x = Age, y = AME, colour = Treatment)) + 
  geom_line(linetype = "solid", size = .8) +
  geom_ribbon(aes(ymin = lower95, ymax = upper95, fill = Treatment), alpha = 0.2) +
  facet_wrap(~Treatment, shrink = TRUE)+
  theme_classic() + 
  scale_color_manual(values = c("Simplification" = "dimgrey", 
                                "Simplification + Social Norm" = "dimgrey",
                                "Simplification + Deterrence" = "dimgrey", 
                                "Simplification + Public Good Appeal" = "dimgrey" 
                              )) + 
  scale_fill_manual(values = c("Simplification" = "dimgrey", 
                               "Simplification + Social Norm" = "dimgrey",
                               "Simplification + Deterrence" = "dimgrey", 
                               "Simplification + Public Good Appeal" = "dimgrey" 
                               )) + 
  labs(y = expression(atop("Predicted Treatment Effect (95% CI)", italic("Property Tax Payment (single payment or instalments)"), 
      )),  x = "Age") +
  geom_hline(yintercept = 0, color = "darkgrey", linetype = "dashed") +
  labs(colour = "Treatment", fill = "Treatment") + 
  xlim(c(20, 100))+  
  theme(legend.position = "none")

print(plot_inter_age_cont)
# ggsave(plot_inter_age_cont, width = 20, height = 12, units = c("cm"), file ="plot_inter_age_cont.jpg")
```

```{r}
## Figure C1c: Heterogeneous Treatment Effects by Age (linear)
model_HTE_age2 <- lm_robust(Tax_paid_install ~ Treatment_full * Age + Tax_total + Tax_declaration_21 + Gender, 
                           data = data_subset_filtered, 
                           se_type = "HC3")

model_HTE_age_inc2 <- margins(model_HTE_age2, 
                         variables = "Treatment_full", 
                         at = list(Age = seq(min(data_subset_filtered$Age, na.rm = TRUE), 
                                             max(data_subset_filtered$Age, na.rm = TRUE), 
                                             by = 1)))

model_HTE_age_inc2 <- as.data.frame(summary(model_HTE_age_inc2))

model_HTE_age_inc2$upper95 <- model_HTE_age_inc2$AME + (model_HTE_age_inc2$SE * 1.96)
model_HTE_age_inc2$lower95 <- model_HTE_age_inc2$AME - (model_HTE_age_inc2$SE * 1.96)

model_HTE_age_inc2 <- model_HTE_age_inc2 %>%  
  mutate(Treatment = case_when(
    factor == "Treatment_fullSimplification" ~ "Simplification", 
    factor == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    factor == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    factor == "Treatment_fullSimplification + Social Norm" ~ "Simplification + Social Norm",
    TRUE ~ factor 
  ))

model_HTE_age_inc2$Treatment <- factor(model_HTE_age_inc2$Treatment, levels = c(
                     "Simplification",
                     "Simplification + Social Norm",
                     "Simplification + Deterrence",
                      "Simplification + Public Good Appeal"
                     ))

plot_inter_age_cont2 <- model_HTE_age_inc2 %>% 
  mutate(AME = AME * 100,
         lower95 = lower95 * 100,
         upper95 = upper95 * 100) %>%
  ggplot(aes(x = Age, y = AME, colour = Treatment)) + 
  geom_line(linetype = "solid", size = .8) +
  geom_ribbon(aes(ymin = lower95, ymax = upper95, fill = Treatment), alpha = 0.2) +
  facet_wrap(~Treatment, shrink = TRUE)+
  theme_classic() + 
  scale_color_manual(values = c("Simplification" = "dimgrey", 
                                "Simplification + Social Norm" = "dimgrey",
                                "Simplification + Deterrence" = "dimgrey", 
                                "Simplification + Public Good Appeal" = "dimgrey" 
                              )) + 
  scale_fill_manual(values = c("Simplification" = "dimgrey", 
                               "Simplification + Social Norm" = "dimgrey",
                               "Simplification + Deterrence" = "dimgrey", 
                               "Simplification + Public Good Appeal" = "dimgrey" 
                               )) + 
  labs(y = expression(atop("Predicted Treatment Effect (95% CI)", italic("Property Tax Payment (single payment or instalments)"), 
      )),  x = "Age") +
  geom_hline(yintercept = 0, color = "darkgrey", linetype = "dashed") +
  labs(colour = "Treatment", fill = "Treatment") + 
  xlim(c(20, 100))+  
  theme(legend.position = "none")

print(plot_inter_age_cont2)
```

```{r}
## Figure C.1d: Predicted Probability of Tax Compliance across Treatment Groups and Age 
plot1_age <- plot_model(model_HTE_age, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  theme(legend.position = "bottom")

plot2_age <- plot_model(model_HTE_age, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Social Norm]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  theme(legend.position = "bottom")

plot3_age <- plot_model(model_HTE_age, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Deterrence]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  theme(legend.position = "bottom")

plot4_age <- plot_model(model_HTE_age, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Public Good Appeal]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  theme(legend.position = "bottom")

grid.arrange(plot1_age, plot2_age, plot3_age, plot4_age)
```

```{r}
## Figure C.1e: Predicted Tax Compliance by Treatment, Age and Gender  
data_subset_filtered <- data_subset_filtered %>% 
  mutate(Gender2 = if_else(Gender == "M", "Male", "Female"))

model_HTE_gender_age <- lm_robust(Tax_paid_install ~ Treatment_full * Gender2 * Age + Treatment_full * Gender2 * I(Age^2) + Tax_total + Tax_declaration_21, 
                           data = data_subset_filtered, 
                           se_type = "HC3")

plot_gender_age <- plot_model(model_HTE_gender_age, 
           type = "pred",
           terms = c("Age[all]",  "Gender2[all]", "Treatment_full[Simplification, Simplification + Social Norm, Simplification + Deterrence, Simplification + Public Good Appeal]]"),
          colors = "bw",
           legend.title = "Gender")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  xlim(c(20, 100)) 

print(plot_gender_age)
```

#### Heterogeneous Treatment Effects - Overall Tax Burden
```{r}
## Figure 4. Heterogeneous Treatment Effects by Tax Burden 
tax_Q1 <- quantile(data_subset_filtered$Tax_total, 0.25)
tax_Q2 <- quantile(data_subset_filtered$Tax_total, 0.50)
tax_Q3 <- quantile(data_subset_filtered$Tax_total, 0.75)

data_subset_filtered <- data_subset_filtered %>% 
  mutate(Tax_quartile = if_else(Tax_total <= tax_Q1, "Q1",
                                if_else(Tax_total > tax_Q1 & Tax_total <= tax_Q2, "Q2", if_else(Tax_total > tax_Q2 & Tax_total <= tax_Q3, "Q3", "Q4"))))

data_subset_filtered$Tax_quartile <-  factor(data_subset_filtered$Tax_quartile,
                                             levels = c("Q1", "Q2", "Q3", "Q4"))

model_HTE_tax <- lm_robust(Tax_paid_install ~ Treatment_full * Tax_quartile + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_tax_inc <- margins(model_HTE_tax, 
                         variables = "Treatment_full", 
                         at = list(Tax_quartile = c("Q1", "Q2", "Q3", "Q4")))

model_HTE_tax_inc <- as.data.frame(summary(model_HTE_tax_inc))

model_HTE_tax_inc$Treatment <- model_HTE_tax_inc$factor
model_HTE_tax_inc$upper95 <- model_HTE_tax_inc$AME + (model_HTE_tax_inc$SE*1.96)
model_HTE_tax_inc$lower95 <- model_HTE_tax_inc$AME - (model_HTE_tax_inc$SE*1.96)

model_HTE_tax_inc <- model_HTE_tax_inc %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    )) %>% 
  mutate(Tax_quartile = case_when(
    Tax_quartile == "Q1" ~ "Q1 (3 - 35 EUR)",
    Tax_quartile == "Q2" ~ "Q2 (35 - 63 EUR)",
    Tax_quartile == "Q3" ~ "Q3 (63 - 99 EUR)",
    Tax_quartile == "Q4" ~ "Q4 (99 - 709 EUR)",
    TRUE ~ Tax_quartile
  ))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_tax <- model_HTE_tax_inc %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100) %>%
  ggplot(aes(colour = Tax_quartile)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), 
                    xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Tax_quartile), 
             position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2)) + 
  theme_classic() + 
  scale_color_grey("Tax burden (quartile)") + 
  scale_fill_grey("Tax burden (quartile)") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Tax burden (quartile)") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  labs(x = expression(atop("Predicted Treatment Effect (95% CI)", italic("Property Tax Payment (single payment or instalments)"), 
      )),  y = "") 

print(plot_inter_tax)
# ggsave(plot_inter_tax, width = 20, height = 12, units = c("cm"), file ="plot_inter_tax.jpg")
```

```{r}
## Figure C.1f: Predicted Probability of Tax Compliance across Treatment Groups and Tax Quartiles  
plot_int_tax <- plot_model(model_HTE_tax, 
           type = "pred",
           terms = c("Treatment_full[all]", "Tax_quartile[all]"),
           #colors = "bw",
           legend.title = "Tax Quartile")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("Q1 (3 - 35 EUR)",
                               "Q2 (35 - 63 EUR)",
                               "Q3 (63 - 99 EUR)",
                               "Q4 (99 - 709 EUR)"))+
  theme(legend.position = "bottom")

print(plot_int_tax)
```

#### Heterogeneous Treatment Effects - Tax Declaration 

```{r}
## Figure C.1g: Heterogeneous Treatment Effects by Tax Declaration  
model_HTE_tax_decl <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full * Tax_declaration_21 + Tax_total + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_tax_decl_inc <- margins(model_HTE_tax_decl, 
                         variables = "Treatment_full", 
                         at = list(Tax_declaration_21 = c("0", "1")))

model_HTE_tax_decl_inc <- as.data.frame(summary(model_HTE_tax_decl_inc))

model_HTE_tax_decl_inc$Treatment <- model_HTE_tax_decl_inc$factor
model_HTE_tax_decl_inc$upper95 <- model_HTE_tax_decl_inc$AME + (model_HTE_tax_decl_inc$SE*1.96)
model_HTE_tax_decl_inc$lower95 <- model_HTE_tax_decl_inc$AME - (model_HTE_tax_decl_inc$SE*1.96)

model_HTE_tax_decl_inc <- model_HTE_tax_decl_inc %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    ))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_tax_decl <- model_HTE_tax_decl_inc %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100,
         Tax_declaration_21 = if_else(Tax_declaration_21 == "0", "No", "Yes")) %>%
  ggplot(aes(colour = Tax_declaration_21)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), 
                    xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Tax_declaration_21), 
             position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1.5, 4), breaks = round(seq(-1.5, 4, by = .5), 2)) + 
  theme_classic() + 
  scale_color_grey("Tax declaration 2021") + 
  scale_fill_grey("Tax declaration 2021") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Tax declaration 2021") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  labs(y = "", 
       x = "Predicted Treatment Effect (95% CI)")

print(plot_inter_tax_decl)
```

```{r}
## Figure C.1h : Predicted Probability of Tax Compliance across Treatment Groups and Tax Declaration 
plot_int_tax_decl <- plot_model(model_HTE_tax_decl, 
           type = "pred",
           terms = c("Treatment_full[all]", "Tax_declaration_21"),
           #colors = "bw",
           legend.title = "Tax declaration 2021")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("No",
                               "Yes"))

print(plot_int_tax_decl)
```

### Average Treatment Effects (DV = Tax paid on-time in a single payment)
```{r}
## Figure C.2a: Treatment Effects Relative to the Control Group

# Basic Model with robust SE (DV = Tax_paid_ontime_full)
model1_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full, 
                    data = data_subset_filtered, se_type = "HC3")

# Model with covariates and robust SE (DV = Tax_paid_ontime_full)
model2_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_summary1_full <- tidy(model1_full, conf.int = TRUE)
model_summary1_full$model <- "Model1"

model_summary2_full <- tidy(model2_full, conf.int = TRUE)
model_summary2_full$model <- "Model2"

model_all_full <- rbind(model_summary1_full, model_summary2_full) 

model_all_full <- model_all_full %>%
  mutate(
    Treatment = gsub("Treatment_full", "", term)
  ) %>% 
  relocate(Treatment, .after = term) 

model_all_full <- model_all_full %>% 
  filter(!grepl("Age|Gender|Tax", term)) %>% 
  mutate(Treatment = if_else(Treatment == "(Intercept)", "Control", Treatment)) %>% 
  select(-term)

model_all_full$Treatment <- factor(model_all_full$Treatment, levels = 
                                c("Simplification + Public Good Appeal",
                                   "Simplification + Deterrence", 
                                   "Simplification + Social Norm", 
                                  "Simplification", 
                                  "Control"
                                  ))

plot_ATE_full <- model_all_full %>% 
  filter(!Treatment == "Control") %>% 
  mutate(estimate = estimate*100,
         conf.low = conf.low*100,
         conf.high = conf.high*100) %>% 
  ggplot(aes(y = Treatment, x = estimate, colour = model)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width=0, position= position_dodge(width=0.5)) +
   geom_point(aes(x = estimate, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2))+
  geom_vline(xintercept = 0, linetype = "dashed")+
  labs(x = "Predicted Treatment Effect (95% CI)", y = "") +
  theme_classic() +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates"), start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(plot_ATE_full)
```

```{r}
## Figure C.2b: Predicted Probability of Tax Compliance across the Control and Treatment Groups (95% CI)
pred_model_1_full <- prediction(model1_full, data_subset_filtered)
pred_model_2_full <- prediction(model2_full, data_subset_filtered)

sum_model1_full <- plyr::ddply(pred_model_1_full, c("Treatment_full"), summarise,
                   mean = mean(fitted, na.rm=TRUE),
                   mean_se = mean(se.fitted, na.rm=TRUE))
sum_model1_full$model <- "Model1"

sum_model2_full <- plyr::ddply(pred_model_2_full, c("Treatment_full"), summarise,
                   mean = mean(fitted, na.rm=TRUE),
                   mean_se = mean(se.fitted, na.rm=TRUE))
sum_model2_full$model <- "Model2"

sum_all_full <- rbind(sum_model1_full, sum_model2_full)

plot_pred_prob_full <- ggplot(sum_all_full, aes(Treatment_full, mean, colour = model)) +
  geom_errorbar(aes(ymin=mean - (1.96 * mean_se), 
                    ymax= mean + (1.96 * mean_se)), width=0, position= position_dodge(width=0.5)) + 
  geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_y_continuous(limits=c(.80, .84))+ 
  ylab("Predicted Probability of Tax Compliance (95% CI)") +
  xlab("")+ 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", 
                                                  "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", 
                                                  "Covariates"), 
                    start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") +
  labs(caption= "")

print(plot_pred_prob_full)
```

### Heterogeneous Treatment Effects (DV = Tax paid on-time in full)
#### Heterogeneous Treatment Effects - Gender
```{r}
## Figure C.2c: Heterogeneous Treatment Effects by Gender  
model_HTE_gender_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full * Gender + Tax_total + Tax_declaration_21 + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_gender_inc_full <- margins(model_HTE_gender_full, 
                         variables = "Treatment_full", 
                         at = list(Gender = c("M", "F")))

model_HTE_gender_inc_full <- as.data.frame(summary(model_HTE_gender_inc_full))

model_HTE_gender_inc_full$Treatment <- model_HTE_gender_inc_full$factor
model_HTE_gender_inc_full$upper95 <- model_HTE_gender_inc_full$AME + (model_HTE_gender_inc_full$SE*1.96)
model_HTE_gender_inc_full$lower95 <- model_HTE_gender_inc_full$AME - (model_HTE_gender_inc_full$SE*1.96)

model_HTE_gender_inc_full <- model_HTE_gender_inc_full %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    )) %>% 
  mutate(Gender = case_when(
    Gender == "F" ~ "Female", 
    Gender == "M" ~ "Male",
    TRUE ~ Gender))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_gender_full <- model_HTE_gender_inc_full %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100) %>% 
  ggplot(aes(colour = Gender)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Gender), position = position_dodge(width=0.5), size=2) +
  theme_classic() + 
  scale_color_grey("Gender") + 
  scale_fill_grey("Gender") + 
  ylab("")+ 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Gender") +
  xlab("Predicted Treatment Effect (95% CI)") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2))

print(plot_inter_gender_full)
```

```{r}
## Figure C.2d: Predicted Probability of Tax Compliance across Treatment Groups and Gender
plot_int_gender_full <- plot_model(model_HTE_gender_full, 
           type = "pred",
           terms = c("Treatment_full[all]", "Gender[all]"),
           #colors = "bw",
           legend.title = "Gender")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("Male", "Female"))

print(plot_int_gender_full)
```

#### Heterogeneous Treatment Effects - Age
```{r}
## Figure C.2e: Heterogeneous Treatment Effects by Age 
model_HTE_age_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full * Age + Treatment_full * I(Age^2) + Tax_total + Tax_declaration_21 + Gender, 
                           data = data_subset_filtered, 
                           se_type = "HC3")

model_HTE_age_inc_full <- margins(model_HTE_age_full, 
                         variables = "Treatment_full", 
                         at = list(Age = seq(min(data_subset_filtered$Age, na.rm = TRUE), 
                                             max(data_subset_filtered$Age, na.rm = TRUE), 
                                             by = 1)))

model_HTE_age_inc_full <- as.data.frame(summary(model_HTE_age_inc_full))

model_HTE_age_inc_full$upper95 <- model_HTE_age_inc_full$AME + (model_HTE_age_inc_full$SE * 1.96)
model_HTE_age_inc_full$lower95 <- model_HTE_age_inc_full$AME - (model_HTE_age_inc_full$SE * 1.96)

model_HTE_age_inc_full <- model_HTE_age_inc_full %>%  
  mutate(Treatment = case_when(
    factor == "Treatment_fullSimplification" ~ "Simplification", 
    factor == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    factor == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    factor == "Treatment_fullSimplification + Social Norm" ~ "Simplification + Social Norm",
    TRUE ~ factor 
  ))

model_HTE_age_inc_full$Treatment <- factor(model_HTE_age_inc_full$Treatment, levels = c(
                     "Simplification",
                     "Simplification + Social Norm",
                     "Simplification + Deterrence",
                      "Simplification + Public Good Appeal"
                     ))

plot_inter_age_cont_full <- model_HTE_age_inc_full %>% 
  mutate(AME = AME * 100,
         lower95 = lower95 * 100,
         upper95 = upper95 * 100) %>%
  ggplot(aes(x = Age, y = AME, colour = Treatment)) + 
  geom_line(linetype = "solid", size = .8) +
  geom_ribbon(aes(ymin = lower95, ymax = upper95, fill = Treatment), alpha = 0.2) +
  facet_wrap(~Treatment, shrink = TRUE)+
  theme_classic() + 
  scale_color_manual(values = c("Simplification" = "dimgrey", 
                                "Simplification + Social Norm" = "dimgrey",
                                "Simplification + Deterrence" = "dimgrey", 
                                "Simplification + Public Good Appeal" = "dimgrey" 
                              )) + 
  scale_fill_manual(values = c("Simplification" = "dimgrey", 
                               "Simplification + Social Norm" = "dimgrey",
                               "Simplification + Deterrence" = "dimgrey", 
                               "Simplification + Public Good Appeal" = "dimgrey" 
                               )) + 
  ylab("Predicted Treatment Effect (95% CI)") + 
  xlab("Age") +
  geom_hline(yintercept = 0, color = "darkgrey", linetype = "dashed") +
  labs(colour = "Treatment", fill = "Treatment") + 
  xlim(c(20, 100))+  
  theme(legend.position = "none")

print(plot_inter_age_cont_full)
```

```{r}
## Figure C.2f: Predicted Probability of Tax Compliance across Treatment Groups and Age
plot1_age_full <- plot_model(model_HTE_age_full, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")

plot2_age_full <- plot_model(model_HTE_age_full, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Social Norm]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")

plot3_age_full <- plot_model(model_HTE_age_full, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Deterrence]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")

plot4_age_full <- plot_model(model_HTE_age_full, 
           type = "pred",
           terms = c("Age[all]", "Treatment_full[Control, Simplification + Public Good Appeal]"),
           colors = "bw",
           legend.title = "Treatment type")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")

grid.arrange(plot1_age_full, plot2_age_full, plot3_age_full, plot4_age_full)
```

```{r}
## Figure C.2g: Predicted Tax Compliance by Treatment, Age and Gender  
data_subset_filtered <- data_subset_filtered %>% 
  mutate(Gender2 = if_else(Gender == "M", "Male", "Female"))

model_HTE_gender_age_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full * Gender2 * Age + Treatment_full * Gender2 * I(Age^2) + Tax_total + Tax_declaration_21, 
                           data = data_subset_filtered, 
                           se_type = "HC3")

plot_gender_age_full <- plot_model(model_HTE_gender_age_full, 
           type = "pred",
           terms = c("Age[all]",  "Gender2[all]", "Treatment_full[Simplification, Simplification + Social Norm, Simplification + Deterrence, Simplification + Public Good Appeal]]"),
          colors = "bw",
           legend.title = "Gender")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "")+
  xlim(c(20, 100)) 

print(plot_gender_age_full)
# ggsave(plot_gender_age, width = 20, height = 12, units = c("cm"), file ="plot_gender_age.jpg")
```

#### Heterogeneous Treatment Effects - Overall Tax Burden

```{r}
## Figure C.2h: Heterogeneous Treatment Effects by Tax Quartile  
tax_Q1 <- quantile(data_subset_filtered$Tax_total, 0.25)
tax_Q2 <- quantile(data_subset_filtered$Tax_total, 0.50)
tax_Q3 <- quantile(data_subset_filtered$Tax_total, 0.75)

data_subset_filtered <- data_subset_filtered %>% 
  mutate(Tax_quartile = if_else(Tax_total <= tax_Q1, "Q1",
                                if_else(Tax_total > tax_Q1 & Tax_total <= tax_Q2, "Q2", if_else(Tax_total > tax_Q2 & Tax_total <= tax_Q3, "Q3", "Q4"))))

data_subset_filtered$Tax_quartile <-  factor(data_subset_filtered$Tax_quartile,
                                             levels = c("Q1", "Q2", "Q3", "Q4"))

model_HTE_tax_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full * Tax_quartile + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_tax_inc_full <- margins(model_HTE_tax_full, 
                         variables = "Treatment_full", 
                         at = list(Tax_quartile = c("Q1", "Q2", "Q3", "Q4")))

model_HTE_tax_inc_full <- as.data.frame(summary(model_HTE_tax_inc_full))

model_HTE_tax_inc_full$Treatment <- model_HTE_tax_inc_full$factor
model_HTE_tax_inc_full$upper95 <- model_HTE_tax_inc_full$AME + (model_HTE_tax_inc_full$SE*1.96)
model_HTE_tax_inc_full$lower95 <- model_HTE_tax_inc_full$AME - (model_HTE_tax_inc_full$SE*1.96)

model_HTE_tax_inc_full <- model_HTE_tax_inc_full %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    )) %>% 
  mutate(Tax_quartile = case_when(
    Tax_quartile == "Q1" ~ "Q1 (3 - 35 EUR)",
    Tax_quartile == "Q2" ~ "Q2 (35 - 63 EUR)",
    Tax_quartile == "Q3" ~ "Q3 (63 - 99 EUR)",
    Tax_quartile == "Q4" ~ "Q4 (99 - 709 EUR)",
    TRUE ~ Tax_quartile
  ))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_tax_full <- model_HTE_tax_inc_full %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100) %>%
  ggplot(aes(colour = Tax_quartile)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), 
                    xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Tax_quartile), 
             position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2)) + 
  theme_classic() + 
  scale_color_grey("Tax burden (quartile)") + 
  scale_fill_grey("Tax burden (quartile)") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Tax burden (quartile)") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  labs(y = "", 
       x = "Predicted Treatment Effect (95% CI)")

print(plot_inter_tax_full)

# ggsave(plot_inter_tax, width = 20, height = 12, units = c("cm"), file ="plot_inter_tax.jpg")
```

```{r}
## Figure C.2i: Predicted Probability of Tax Compliance across Treatment Groups and Tax Quartiles 
plot_int_tax_full <- plot_model(model_HTE_tax_full, 
           type = "pred",
           terms = c("Treatment_full[all]", "Tax_quartile[all]"),
           #colors = "bw",
           legend.title = "Tax Quartile")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("Q1 (3 - 35 EUR)",
                               "Q2 (35 - 63 EUR)",
                               "Q3 (63 - 99 EUR)",
                               "Q4 (99 - 709 EUR)"))

print(plot_int_tax_full)
```

#### Heterogeneous Treatment Effects - Tax Declaration 

```{r}
## Figure C.2j: Heterogeneous Treatment Effects by Tax Declaration  
model_HTE_tax_decl_full <- lm_robust(Tax_paid_ontime_full ~ Treatment_full * Tax_declaration_21 + Tax_total + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_HTE_tax_decl_inc_full <- margins(model_HTE_tax_decl_full, 
                         variables = "Treatment_full", 
                         at = list(Tax_declaration_21 = c("0", "1")))

model_HTE_tax_decl_inc_full <- as.data.frame(summary(model_HTE_tax_decl_inc_full))

model_HTE_tax_decl_inc_full$Treatment <- model_HTE_tax_decl_inc_full$factor
model_HTE_tax_decl_inc_full$upper95 <- model_HTE_tax_decl_inc_full$AME + (model_HTE_tax_decl_inc_full$SE*1.96)
model_HTE_tax_decl_inc_full$lower95 <- model_HTE_tax_decl_inc_full$AME - (model_HTE_tax_decl_inc_full$SE*1.96)

model_HTE_tax_decl_inc_full <- model_HTE_tax_decl_inc_full %>% 
  mutate(Treatment = case_when(
    Treatment == "Treatment_fullSimplification" ~ "Simplification", 
    Treatment == "Treatment_fullSimplification + Deterrence" ~ "Simplification + Deterrence", 
    Treatment == "Treatment_fullSimplification + Public Good Appeal" ~ "Simplification + Public Good Appeal", 
    Treatment == "Treatment_fullSimplification + Social Norm"  ~ "Simplification + Social Norm",
    TRUE ~ Treatment 
    ))

treatment_order <- c(
                     "Simplification + Public Good Appeal",
                     "Simplification + Deterrence",
                     "Simplification + Social Norm",
                     "Simplification"
                     )

plot_inter_tax_decl_full <- model_HTE_tax_decl_inc_full %>% 
  mutate(AME = AME*100,
         lower95 = lower95*100,
         upper95 = upper95*100,
         Tax_declaration_21 = if_else(Tax_declaration_21 == "0", "No", "Yes")) %>%
  ggplot(aes(colour = Tax_declaration_21)) + 
  geom_errorbar(aes(y = factor(Treatment, level = treatment_order), 
                    xmin=lower95, xmax=upper95), 
                width=0, position=position_dodge(width=0.5)) + 
  geom_point(aes(y = factor(Treatment, level = treatment_order), x=AME, shape = Tax_declaration_21), 
             position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1.5, 4), breaks = round(seq(-1.5, 4, by = .5), 2)) + 
  theme_classic() + 
  scale_color_grey("Tax declaration 2021") + 
  scale_fill_grey("Tax declaration 2021") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
  scale_shape_discrete("Tax declaration 2021") +
  geom_vline(xintercept=0, color="grey", linetype = "dashed")+
  labs(y = "", 
       x = "Predicted Treatment Effect (95% CI)")

print(plot_inter_tax_decl_full)
```

```{r}
## Figure C.2k: Predicted Probability of Tax Compliance across Treatment Groups and Tax Declaration  
plot_int_tax_decl_full <- plot_model(model_HTE_tax_decl_full, 
           type = "pred",
           terms = c("Treatment_full[all]", "Tax_declaration_21"),
           #colors = "bw",
           legend.title = "Tax declaration 2021")+
  theme_classic() + 
  labs(y = "Predicted Probability of Tax Compliance (95% CI)", title = "", x = "")+ 
  scale_colour_grey(labels = c("No",
                               "Yes"))

print(plot_int_tax_decl_full)
```

### C.3 Full Regression Tables
```{r}
## Table C.3a: Full Regression Table for Figure 1
modelsummary(models = list("Model 1" = model1, "Model 2" = model2),
                           statistic = "std.error",
                           stars = c('***' = .001, '**' = .01, '*' = .05),
             coef_omit = c(10), 
                           coef_rename = c("Control", 
                                        "Simplification", 
                                        "Simplification + Social Norm", 
                                        "Simplification + Deterrence", 
                                        "Simplification + Public Good Appeal",
                                        "Tax total", 
                                        "Tax declaration 2021",
                                        "Female", 
                                        "Age"),
                           output = "kableExtra" 
                           ) %>% 
    kable_styling(full_width = FALSE, font_size = 10, latex_options = "HOLD_position")
```

```{r}
## Table C.3b: Full Regression Table for Figure C.2a
modelsummary(models = list("Model 1" = model1_full, "Model 2" = model2_full),
                           statistic = "std.error",
                           stars = c('***' = .001, '**' = .01, '*' = .05),
                           coef_omit = c(10), 
                           coef_rename = c("Control", 
                                        "Simplification", 
                                        "Simplification + Social Norm", 
                                        "Simplification + Deterrence", 
                                        "Simplification + Public Good Appeal",
                                        "Tax total", 
                                        "Tax declaration 2021",
                                        "Female", 
                                        "Age"),
                           output = "kableExtra") %>% 
    kable_styling(full_width = FALSE, font_size = 10, latex_options = "HOLD_position")
```

### Figure C.4a:** Treatment Effects Relative to the Control Group (Different Heteroskedasticity-Consistent Covariance Matrix estimators)
```{r}
## Figure C.4a:** Treatment Effects Relative to the Control Group (Different Heteroskedasticity-Consistent Covariance Matrix estimators)

# Model with covariates and standard SE (Non-robust)
model_classic <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "classical")

# Model with covariates and robust SE (HC0)
model_HC0 <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC0")

# Model with covariates and robust SE (HC1 - Stata)
model_HC1 <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC1")

# Model with covariates and robust SE (HC2)
model_HC2 <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC2")

# Model with covariates and robust SE (HC3)
model_HC3 <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered, se_type = "HC3")

model_summary_classic <- tidy(model_classic, conf.int = TRUE)
model_summary_classic$model <- "Classic"

model_summary_HC0 <- tidy(model_HC0, conf.int = TRUE)
model_summary_HC0$model <- "HC0"

model_summary_HC1 <- tidy(model_HC1, conf.int = TRUE)
model_summary_HC1$model <- "HC1"

model_summary_HC2 <- tidy(model_HC2, conf.int = TRUE)
model_summary_HC2$model <- "HC2"

model_summary_HC3 <- tidy(model_HC3, conf.int = TRUE)
model_summary_HC3$model <- "HC3"


model_all_full <- rbind(model_summary_classic, model_summary_HC0, model_summary_HC1, model_summary_HC2, model_summary_HC3) 

model_all_full <- model_all_full %>%
  mutate(
    Treatment = gsub("Treatment_full", "", term)
  ) %>% 
  relocate(Treatment, .after = term) 

model_all_full <- model_all_full %>% 
  filter(!grepl("Age|Gender|Tax", term)) %>% 
  mutate(Treatment = if_else(Treatment == "(Intercept)", "Control", Treatment)) %>% 
  select(-term)

model_all_full$Treatment <- factor(model_all_full$Treatment, levels = 
                                c("Simplification + Public Good Appeal",
                                   "Simplification + Deterrence", 
                                   "Simplification + Social Norm", 
                                  "Simplification", 
                                  "Control"
                                  ))

plot_SE_robust <- model_all_full %>% 
  filter(!Treatment == "Control") %>% 
  mutate(estimate = estimate*100,
         conf.low = conf.low*100,
         conf.high = conf.high*100) %>% 
  ggplot(aes(y = Treatment, x = estimate, colour = model)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width=0, position= position_dodge(width=0.5)) +
   geom_point(aes(x = estimate, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2))+
  geom_vline(xintercept = 0, linetype = "dashed")+
  labs(x = "Predicted Treatment Effect (95% CI)", y = "") +
  theme_classic() +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("Non-robust SE", "HC0 estimator", "HC1 - Stata estimator", "HC2 estimator", "HC3 estimator")) +
  scale_colour_grey(name = "Models", labels = c("Non-robust SE", "HC0 estimator", "HC1 - Stata estimator", "HC2 estimator", "HC3 estimator"), start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 
```

### D Regression Results under the Alternative Outlier Removel Process (2.5 percentile rule) and Full Dataset

#### D.1 Regression Results under the Alternative Outlier Removel Process (2.5 percentile rule)
```{r}
# 2.5 and 97.5 percentile bounds 
lower_bound <- quantile(data_subset$Tax_total, 0.025)
upper_bound <- quantile(data_subset$Tax_total, 0.975)

data_subset <- data_subset %>% 
  mutate(perc_out = if_else(Tax_total > upper_bound | Tax_total < lower_bound, 1, 0))
  
data_subset_filtered2 <- data_subset %>% 
  filter(perc_out != 1)

## Figure D.1a 
# Basic Model with robust SE
model1b <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full, 
                    data = data_subset_filtered2, se_type = "HC3")

# Model with covariates and robust SE
model2b <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset_filtered2, se_type = "HC3")

model_summary1b <- tidy(model1b, conf.int = TRUE)
model_summary1b$model <- "Model1b"

model_summary2b <- tidy(model2b, conf.int = TRUE)
model_summary2b$model <- "Model2b"

model_allb <- rbind(model_summary1b, model_summary2b) 

model_allb <- model_allb %>%
  mutate(
    Treatment = gsub("Treatment_full", "", term)
  ) %>% 
  relocate(Treatment, .after = term) 

model_allb <- model_allb %>% 
  filter(!grepl("Admin|Age|Gender|Tax", term)) %>% 
  mutate(Treatment = if_else(Treatment == "(Intercept)", "Control", Treatment)) %>% 
  select(-term)

model_allb$Treatment <- factor(model_allb$Treatment, levels = 
                                c("Simplification + Public Good Appeal",
                                   "Simplification + Deterrence", 
                                   "Simplification + Social Norm", 
                                  "Simplification", 
                                  "Control"
                                  ))

figure1b <- model_allb %>% 
  filter(!Treatment == "Control") %>% 
  mutate(estimate = estimate*100,
         conf.low = conf.low*100,
         conf.high = conf.high*100) %>% 
  ggplot(aes(y = Treatment, x = estimate, colour = model)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width=0, position= position_dodge(width=0.5)) +
   geom_point(aes(x = estimate, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2))+
  geom_vline(xintercept = 0, linetype = "dashed")+
  labs(x = "Predicted Treatment Effect (95% CI)", y = "") +
  theme_classic() +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates"), start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(figure1b)
```

```{r}
## Figure D.1b: Predicted Probability of Tax Compliance across the Control and Treatment Groups (95% CI)
pred_model_1b <- prediction(model1b, data_subset_filtered)
pred_model_2b <- prediction(model2b, data_subset_filtered)

sum_model1b <- plyr::ddply(pred_model_1b, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model1b$model <- "Model1b"

sum_model2b <- plyr::ddply(pred_model_2b, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model2b$model <- "Model2b"

sum_allb <- rbind(sum_model1b, sum_model2b)

plot_pred_probb <- ggplot(sum_allb, aes(Treatment_full, mean, colour = model)) +
  geom_errorbar(aes(ymin=mean - (1.96 * mean_se), 
                    ymax= mean + (1.96 * mean_se)), width=0, position= position_dodge(width=0.5)) + 
  geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_y_continuous(limits=c(.92, .95))+ 
  ylab("Predicted Probability of Tax Compliance (95% CI)") +
  xlab("")+ 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", 
                                                  "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", 
                                                  "Covariates"), 
                    start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(plot_pred_probb)
```

#### D.2 Regression Results with the Untrimmed Dataset
```{r}
## Figure D.2a 
# Basic Model with robust SE
model1c <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full, 
                    data = data_subset, se_type = "HC3")

# Model with covariates and robust SE
model2c <- lm_robust(Tax_paid_ontime_full_or_part ~ Treatment_full + Tax_total + Tax_declaration_21 + Gender + Age + I(Age^2), data = data_subset, se_type = "HC3")

model_summary1c <- tidy(model1c, conf.int = TRUE)
model_summary1c$model <- "Model1c"

model_summary2c <- tidy(model2c, conf.int = TRUE)
model_summary2c$model <- "Model2c"

model_allc <- rbind(model_summary1c, model_summary2c) 

model_allc <- model_allc %>%
  mutate(
    Treatment = gsub("Treatment_full", "", term)
  ) %>% 
  relocate(Treatment, .after = term) 

model_allc <- model_allc %>% 
  filter(!grepl("Admin|Age|Gender|Tax", term)) %>% 
  mutate(Treatment = if_else(Treatment == "(Intercept)", "Control", Treatment)) %>% 
  select(-term)

model_allc$Treatment <- factor(model_allc$Treatment, levels = 
                                c("Simplification + Public Good Appeal",
                                   "Simplification + Deterrence", 
                                   "Simplification + Social Norm", 
                                  "Simplification", 
                                  "Control"
                                  ))

figure1c <- model_allc %>% 
  filter(!Treatment == "Control") %>% 
  mutate(estimate = estimate*100,
         conf.low = conf.low*100,
         conf.high = conf.high*100) %>% 
  ggplot(aes(y = Treatment, x = estimate, colour = model)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width=0, position= position_dodge(width=0.5)) +
   geom_point(aes(x = estimate, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_x_continuous(limits=c(-1, 4), breaks = round(seq(-1, 4, by = .5), 2))+
  geom_vline(xintercept = 0, linetype = "dashed")+
  labs(x = "Predicted Treatment Effect (95% CI)", y = "") +
  theme_classic() +
  theme(axis.text.y = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates"), start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(figure1c)
```

```{r}
## Figure D.2b: Predicted Probability of Tax Compliance across the Control and Treatment Groups (95% CI)
pred_model_1c <- prediction(model1c, data_subset_filtered)
pred_model_2c <- prediction(model2c, data_subset_filtered)

sum_model1c <- plyr::ddply(pred_model_1c, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model1c$model <- "Model1c"

sum_model2c <- plyr::ddply(pred_model_2c, 
                          c("Treatment_full"), 
                          summarise,
                          mean = mean(fitted, na.rm=TRUE),
                          mean_se = mean(se.fitted, na.rm=TRUE))
sum_model2c$model <- "Model2c"

sum_allc <- rbind(sum_model1c, sum_model2c)

plot_pred_probc <- ggplot(sum_allc, aes(Treatment_full, mean, colour = model)) +
  geom_errorbar(aes(ymin=mean - (1.96 * mean_se), 
                    ymax= mean + (1.96 * mean_se)), width=0, position= position_dodge(width=0.5)) + 
  geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
  scale_y_continuous(limits=c(.92, .95))+ 
  ylab("Predicted Probability of Tax Compliance (95% CI)") +
  xlab("")+ 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8)) + 
 scale_shape_discrete(name = "Models", labels = c("No covariates", 
                                                  "Covariates")) +
  scale_colour_grey(name = "Models", labels = c("No covariates", 
                                                  "Covariates"), 
                    start = 0.2, end = 0.8, aesthetics = "colour") +
  theme(legend.position="right") 

print(plot_pred_probc)
```