Sample Codes

[shoonlee]

Texreg

texreg(list(reg_perish_kg_stack, reg_store_kg_stack, reg_freq, reg_dining), float.pos = '!t', caption.above = T,
       caption = "Margins of Grocery Purchase Adjustment", digits = 3, leading.zero = F,
       include.adjrs = F, include.rs = F, include.ci = F, include.rmse = F,
       custom.model.names = c("(1)", "(2)", "(3)", "(4)"),
       custom.gof.rows = list('Dep.Var' = c("Perishable KG", "Storable KG", "Shopping Freq", "Eatout Freq")),
       stars = NULL,
       custom.coef.map = list('post:ups_group' = 'Post $\\times$ Treated', 'p_waste3' = 'Fee (Cents)'),
       label = "table:food_margin", fontsize = "small", use.packages = FALSE, include.groups = FALSE,
       custom.note = "\\parbox{.75\\linewidth}
{\\vspace{2pt} Note: This table is produced from equation (\\ref{eq:food}). Dependent variables in columns (1) and (2) are the total annual grocery purchase amount per household. In columns (3) and (4), annual perishable and storable grocery purchase amount per household are regressed, respectively. UPS is a dummy variable that takes 1 if a household is subject to the UPS. Fee (Cents) is per KG waste fee converted to the US Cents. All standard errors are clustered at the municipality level.\\\\}")

[shoonlee]

Creating a plot (non-parametric damage function). You can ignore the back-transformation part and start directly from tidying regression results.

# Backtransformation using Delta method (irrelevant)

delta_exp3 <- function(case, i, model){

  if(case == "ctrl_pre"){
    tmp1 <- paste0("fl_", i)
    tmp_combined <- paste0("exp(0.686/2)*exp(", tmp1, ")")
  } else if(case == "ctrl_post"){
    tmp1 <- paste0("fl_", i)
    tmp2 <- paste0("postfl_", i)
    tmp_combined <- paste0("exp(0.686/2)*exp(", paste(c(tmp1, tmp2), collapse = "+"), ")")
  } else if(case == "disc_pre"){
    tmp1 <- paste0("fl_", i)
    tmp2 <- paste0("discfl_", i)
    tmp_combined <- paste0("exp(0.686/2)*exp(", paste(c(tmp1, tmp2), collapse = "+"), ")")
  } else{
    tmp1 <- paste0("fl_", i)
    tmp2 <- paste0("discfl_", i)
    tmp3 <- paste0("postfl_", i)
    tmp4 <- paste0("pdfl_", i)
    tmp_combined <- paste0("exp(0.686/2)*exp(", paste(c(tmp1, tmp2, tmp3, tmp4), collapse = "+"), ")")
  }

  deltaMethod(get(model), tmp_combined)    

}

case <- rep(rep(c("ctrl_pre", "ctrl_post", "disc_pre", "disc_post"), each = 5), 5)
i <- rep(c(2:6), 20)
model <- rep(c("regpop_yr_comm", "regpop_yr_comm_high", 
               "regpop_yr_comm_low", "regavgdmg_yr_comm", "regct_yr_comm"), each = 20)

delta_tmp <- pmap(.l = list(case, i, model), .f = delta_exp3)

# Tidy the regression results
delta_df_all <- rbindlist(delta_tmp) %>%
  rename(conf.low = `2.5 %`, conf.high = `97.5 %`, estimate = Estimate) %>%
  mutate(flsize = rep(c("2-10", "10-20", "20-30", "30-40", "40-50"), 20),
         term = case,
         model = model) %>%
  select(-SE) 

# Make model choice
delta_df <- delta_df_all %>%
  filter(model == "regpop_yr_comm")

  desired_order <- c("2-10", "10-20", "20-30", "30-40", "40-50")

  delta_df$flsize <- factor( as.character(delta_df$flsize), levels=desired_order )

  delta_df_ctrl <- delta_df %>% slice(1:10)
  delta_df_ctrl$term <- factor(delta_df_ctrl$term, levels=c("ctrl_pre", "ctrl_post"), labels = c("Pre", "Post"))
  delta_df_disc <- delta_df %>% slice(11:20)
  delta_df_disc$term <- factor(delta_df_disc$term, levels=c("disc_pre", "disc_post"), labels = c("Pre", "Post"))

  # w/o SE

  ggplot(data = delta_df_ctrl, aes(x=flsize, y=estimate, group = term, lty=term, shape=term)) + 
    geom_line() + 
    geom_point() + 
    theme_minimal() + xlab("Flood Size") + ylab("Per Capita Damage to the Baseline") +
    theme(text = element_text(size=13), panel.grid.minor = element_blank(), legend.position="bottom") + 
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0)