# '
# Creates the summary statistics for health outcome and control variables
# as Supplementary Table SI2 of the paper.
# Uses panel data stored to "bld/data",
# variables names from "src/code/specifications.json"
# and saves the plot as a tex file to "bld/tables/".
# 
# '

rm(list = ls())
options(scipen = 999) # Turn off scientific notation
options(warn = 0) # 0 to turn on again, -1 to turn off

# Load libraries
library(tidyverse) # More intuitive data wrangling
library(xtable) # Create flexible tex files
library(rjson) # Read JSON files in R

# Load model specification
model <- fromJSON(file = paste0("src/code/specifications.json"))

# Load analysis data
analysis <- readRDS(paste0("bld/data/analysis_year.RDS"))

# Keep only obs at 100km bandwidth
dat <- subset(analysis, analysis$dist2cutoff <= model$distance_to_border & analysis$year > 2002)

# Write table header
final_table <- data.frame((matrix(nrow = 1, ncol = 5)))
final_table[1] <- c("")
final_table[2] <- c("Outside")
final_table[3] <- c("Inside")
final_table[4] <- c("Outside")
final_table[5] <- c("Inside")


# Controls variables

# Drop duplicates
dat_distinct <- dat %>%
  filter(icd == "J") 

# Before period
d <- subset(dat_distinct, dat_distinct$year < 2007) %>% distinct(mcode, year, .keep_all = TRUE)
# Fill data into table conditional on model
controls_table = data.frame(matrix(nrow = (length(model$controls) - 1), ncol = 5))
controls_table[1] <- sapply(2:length(model$controls), function(i) paste(sub("\\\\ ", "", model$controls[[i]], fixed=TRUE)))
controls_table[2] <- sapply(2:length(model$controls), function(i) paste0(
  round(mean(d[[names(model$controls[i])]][d$treated == 0], na.rm = TRUE), 2),
  " (",
  round(sd(d[[names(model$controls[i])]][d$treated == 0], na.rm = TRUE), 2),
  ")"
))
controls_table[3] <- sapply(2:length(model$controls), function(i) paste0(
  round(mean(d[[names(model$controls[i])]][d$treated == 1], na.rm = TRUE), 2),
  " (",
  round(sd(d[[names(model$controls[i])]][d$treated == 1], na.rm = TRUE), 2),
  ")"
))

# After period
d <- subset(dat_distinct, dat_distinct$year > 2006)
# Fill data into table conditional on model
controls_table[4] <- sapply(2:length(model$controls), function(i) paste0(
  round(mean(d[[names(model$controls[i])]][d$treated == 0], na.rm = TRUE), 2),
  " (",
  round(sd(d[[names(model$controls[i])]][d$treated == 0], na.rm = TRUE), 2),
  ")"
))
controls_table[5] <- sapply(2:length(model$controls), function(i) paste0(
  round(mean(d[[names(model$controls[i])]][d$treated == 1], na.rm = TRUE), 2),
  " (",
  round(sd(d[[names(model$controls[i])]][d$treated == 1], na.rm = TRUE), 2),
  ")"
))

# Health outcome variables

# Before period
d <- subset(dat, dat$year < 2007)
# Fill data into table conditional on model
cases_table = data.frame(matrix(nrow = (length(model$outcomes) - 1), ncol = 5))
cases_table[1] <- sapply(2:length(model$outcomes), function(i) paste(sub("\\\\ ", "", model$outcomes[[i]], fixed=TRUE)))
cases_table[2] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["cases_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["cases_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))
cases_table[3] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["cases_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["cases_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))

# After period
d <- subset(dat, dat$year > 2006)
# Fill data into table conditional on model
cases_table[4] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["cases_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["cases_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))
cases_table[5] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["cases_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["cases_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))

# Before period
d <- subset(dat, dat$year < 2007)
# Fill data into table conditional on model
deaths_table = data.frame(matrix(nrow = (length(model$outcomes) - 1), ncol = 5))
deaths_table[1] <- sapply(2:length(model$outcomes), function(i) paste(sub("\\\\ ", "", model$outcomes[[i]], fixed=TRUE)))
deaths_table[2] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["death_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["death_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))
deaths_table[3] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["death_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["death_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))

# After period
d <- subset(dat, dat$year > 2006)
# Fill data into table conditional on model
deaths_table[4] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["death_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["death_prevalence"]][d$treated == 0 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))
deaths_table[5] <- sapply(2:length(model$outcomes), function(i) paste0(
  round(mean(d[["death_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  " (",
  round(sd(d[["death_prevalence"]][d$treated == 1 & d$icd == names(model$outcomes[i])], na.rm = TRUE), 2),
  ")"
))

# Add headers to each section of the table
control_header = c("\\multicolumn{1}{l}{\\textbf{Control variables}}", "", "", "", "")
cases_header = c("\\multicolumn{1}{l}{\\textbf{Hospitalizations per 1,000}}", "", "", "", "")
deaths_header = c("\\multicolumn{1}{l}{\\textbf{Deaths per 1,000}}", "", "", "", "")

# Add units description to table
controls_table[ ,1] <- gsub('capita', 'capita (1,000R\\\\$/pop)', controls_table[ ,1])
# controls_table[ ,1] <- gsub('share', 'share (\\\\%)', controls_table[ ,1])
controls_table[ ,1] <- gsub('density', 'density (pop/\\$km^2)\\$', controls_table[ ,1])
controls_table[ ,1] <- gsub('Doctors', 'Doctors (per 1,000)', controls_table[ ,1])
controls_table[ ,1] <- gsub('establishments', 'establishments (per 1,000)', controls_table[ ,1])
controls_table[ ,1] <- gsub('Temperature', 'Temperature (\\\\degree C)', controls_table[ ,1])
controls_table[ ,1] <- gsub('Precipitation', 'Precipitation (mm/day)', controls_table[ ,1])
controls_table[ ,1] <- gsub('Humidity', 'Humidity (\\\\%)', controls_table[ ,1])
controls_table[ ,1] <- gsub('speed', 'speed (m/s)', controls_table[ ,1])
controls_table[ ,1] <- gsub('Fines', 'Fines (R\\\\$/pop)', controls_table[ ,1])
controls_table[ ,1] <- gsub('area', 'area (\\\\%)', controls_table[ ,1])
controls_table[ ,1] <- gsub('direction', 'direction (in \\degree)', controls_table[ ,1])

# Append tables
final_table = rbind(final_table, control_header, controls_table, cases_header, cases_table, deaths_header, deaths_table)

# Create latex table
tex_final_table = xtable(
  final_table,
  #caption = "Mean and Standard Deviations of Pre- and Post-ASM Variables by Treatment",
  #label = "table:summary",
  #align = "lrrrrr"
)

# Add note and top row
comment          <- list()
comment$pos      <- as.list(c(0, nrow(final_table)))
comment$command  <- 
  as.vector(
    c(
      paste(
        "\\begin{tabular}{rrrrr} \n",
        "\\\\[-1.8ex]\\hline \n", "\\hline \\\\[-1.8ex] \n",
        "& \\multicolumn{2}{c}{Before} & \\multicolumn{2}{c}{After} \\\\ \n",
        sep = ""
      ),
      paste(
        "\\hline \n", "\\hline \\\\[-1.8ex] \n",
        "\\end{tabular} \n",
        "\\begin{tablenotes}[para,flushleft] \n",
        "\\textit{Note:} Numbers represent averages of each group in each period. The outside group refers to municipalities whose capitals
        lie outside the Amazon biome but no further than 100 km away from the border (correspondingly for the inside group).
        The before period is from 2003 to 2006, and the after period from 2007 to 2017.
        Temperature is in degree Celsius; precipitation is in millimeters per day; 
        humidity is in percent; wind speed in meters per second;
        GDP per capita is in thousand Brazilian Reais in 2012 constant prices;
        agricultural GDP share is the GDP share of the agricultural sector;
        population density is inhabitants per square km; doctors and health establishments are per 1,000 inhabitants;
        and wind direction is in degrees.
        Hospitalizations and deaths are per 1,000 inhabitants associated with the named health problem/disease
        according to the ICD-10 classification. STDs are sexually transmitted diseases. Standard deviations in parentheses. \\\\ \n",
        "\\end{tablenotes} \n",
        sep = ""
      )
    ), mode="character"
)

# Save the latex table
print(
  tex_final_table,
  sanitize.text.function = function(x){x},
  file = paste("bld/tables/summary.tex", sep = "/"),
  include.rownames = FALSE,
  include.colnames = FALSE,
  caption.placement = "top",
  add.to.row = comment,
  hline.after = c(1),
  floating = FALSE,
  only.contents = TRUE
)