# Not run! 
fit <- lm(Y ~ X, data = dat)

tab <- huxreg(fit,
              # Rename coefficients
              coefs = c("Cool variable name" = "X")
              # Add the desired summary statistics
              statistics = c(N = "nobs", R2 = "r.squared"),
              # Add a note
              note = "Note. Important note about this table") %>% 
  # Convert the huxtable table to flextable
  as_flextable()

# Set a custom width to ensure it fits on the page
tab %>% flextable::width(width = 1)

library(tidyverse)

library(dbplyr)

path <- "/Volumes/the safe box/CCD_2001_2016/ccd_2001-2016.sqlite"

ccd <- DBI::dbConnect(RSQLite::SQLite(), path)

DBI::dbListTables(ccd)

query <- 
  "SELECT * 
    FROM (SELECT fips, year AS e_year, leaid AS e_leaid, race, enrollment 
            FROM ccd_district_enroll 
            WHERE (race IN ('White', 'Black', 'Hispanic', 'Asian', 'American Indian or Alaska Native', 'Total') AND year == 2015)) AS e 
    INNER JOIN 
      (SELECT leaid, year, rev_total, rev_fed_total, rev_local_total, rev_state_total 
        FROM ccd_district_finance) AS f 
    ON (e.e_leaid = f.leaid AND e.e_year = f.year)"

dat <- tbl(ccd, sql(query)) %>% 
  collect() %>% 
  select(-starts_with("e_")) %>% 
  pivot_wider(names_from = race, values_from = enrollment) %>% 
  janitor::clean_names() %>% 
  mutate_at(vars(white, black, hispanic, asian), list(per = ~ round(100 * (./total), 2))) %>% 
  na_if(-2) %>% 
  na_if(-1) %>% 
  filter_at(vars(contains("rev")), all_vars(. != 0))

enroll <- get_education_data(level = "school-districts",
                        source = "ccd",
                        topic = "enrollment",
                        filters = list(year = 2015, grade = 99),
                        by = list("race"),
                        add_labels = TRUE)

fin <- get_education_data(level = "school-districts",
                     source = "ccd",
                     topic = "finance",
                     filters = list(year = 2015),
                     add_labels = TRUE)

dat <- enroll %>% 
  select(leaid, year, race, enrollment)
  inner_join(fin %>% 
               select(leaid, year, rev_total, rev_fed_total, rev_local_total),
             by = c("leaid", "year")) %>% 
    pivot_wider(names_from = race, values_from = enrollment) %>% 
    janitor::clean_names()

library(huxtable)
fit1 <- lm(log(rev_local_total) ~ white_per, data = dat)

huxreg(fit1)

huxreg(fit1,
       statistics = c(N = "nobs", R2 = "r.squared"))

huxreg(fit1,
       statistics = c(`Number of observations` = "nobs", `R-squared` = "r.squared"))

huxreg(fit1,
       coefs = c("%White Students" = "white_per"),
       statistics = c(`Number of observations` = "nobs", `R-squared` = "r.squared"))

fit2 <- update(fit1, . ~ . + log(rev_state_total))
fit3 <- update(fit1, . ~ . + log(rev_fed_total))
fit4 <- update(fit2, . ~ . + log(rev_fed_total))

fits <- list(fit1, fit2, fit3, fit4)

huxreg(fits,
       coefs = c("%White Students" = "white_per",
                 "State revenue (logged)" = "log(rev_state_total)",
                 "Federal revenue (logged)" = "log(rev_fed_total)"),
       statistics = c(N = "nobs", R2 = "r.squared"))

fits <- list("Model 1" = fit1, "Model 2" = fit2, 
             "Model 3" = fit3, "Model 4" = fit4)

huxreg(fits,
       statistics = c(N = "nobs", R2 = "r.squared"),
       note = "Note. Data from the Common Core of Data, National Center for Education Statistics")

state_fits <- map(fits, ~ update(.x, . ~ . + fips))

huxreg(state_fits, 
       statistics = c(N = "nobs", R2 = "r.squared"),
       note = "Note. Standard errors clustered at the state level. Data from the Common Core of Data, National Center for Education Statistics")

fips_coefs <- names(coef(state_fits[[1]]))[str_detect(names(coef(state_fits[[1]])), "fips")]

huxreg(state_fits, 
       omit_coefs = fips_coefs,
       statistics = c(N = "nobs", R2 = "r.squared"),
       note = "Note. Standard errors clustered at the state level. Data from the Common Core of Data, National Center for Education Statistics")

huxreg(state_fits, 
       coefs = c("%White Students" = "white_per",
                 "State revenue (logged)" = "log(rev_state_total)",
                 "Federal revenue (logged)" = "log(rev_fed_total)"),
       statistics = c(N = "nobs", R2 = "r.squared"),
       note = "Note. Standard errors clustered at the state level. All models include a control for state. Data from the Common Core of Data, National Center for Education Statistics")

tab <- huxreg(c(fits, state_fits),
         coefs = c("%White Students" = "white_per",
                 "State revenue (logged)" = "log(rev_state_total)",
                 "Federal revenue (logged)" = "log(rev_fed_total)"),
       statistics = c(N = "nobs", R2 = "r.squared"),
       note = "Note. Standard errors clustered at the state level. Data from the Common Core of Data, National Center for Education Statistics")

dim(tab)

tab %>% insert_row(c("State Control", rep("N", 4), rep("Y", 4)), after = 7)

tab %>% 
  insert_row(c("State control", rep("N", 4), rep("Y", 4)), after = 7) %>% 
  set_bottom_border(row = 8, col = 2:9, value = 0)

models <- list(m1, m2, m3, m4, m5, m6)

did_tidy <- map(models, tidy)

did_tidy <-
  did_tidy %>%
  map(~ mutate(.,
               estimate =
                 case_when(
                 p.value <= 0.001 ~ paste0(round(estimate, 3), " ***"),
                 dplyr::between(p.value, 0.001, 0.01) ~ paste0(round(estimate, 3), " **"),
                 dplyr::between(p.value, 0.01, 0.05) ~ paste0(round(estimate, 3), " *"),
                 TRUE ~ as.character(round(estimate, 3))
               ),
               std.error = paste0("(", round(std.error, 3), ")")
               )) %>%
  map(~ filter(.,
               str_detect(term, ":"))
      ) %>%
  map(~ select(.,
               term, estimate, std.error))

did_tidy <-
  did_tidy %>%
  map(., ~
        reshape2::melt(.x,
                       id.var = "term")
      )

did_tidy <-
  reduce(did_tidy,
         full_join,
         by = c("term",
                "variable")
         ) %>%
  select(-variable)

did_tidy[["term"]] <-
  c("Increase %Black or Latinx, 2013-2015 X post-2013",
    "")

tab <-
  did_tidy %>%
  hux() %>%
  insert_row("",
             "Suspension Rate", "", "",
             "Attendance Rate", "", "",
             after = 0) %>%
  set_colspan(1, c(2, 5), 3) %>%
  insert_row("",
             "All grades", "Grades 3-5", "Grades 6-8",
             "All grades", "Grades 3-5", "Grades 6-8",
             after = 1) %>%
  # 4 total rows, 5 total columns here
  insert_row("School Fixed Effects",
             "Y", "Y",
             "Y", "Y",
             "Y", "Y",
             after = 4) %>%
  insert_row("Year Fixed Effects",
             "Y", "Y",
             "Y", "Y",
             "Y", "Y",
             after = 5) %>%
  insert_row("Observations",
             map(models, `[`, "N") %>%
               unlist(),
             after = 6) %>%
  insert_row("R2",
             map(models, glance) %>%
               map(`[`, "r.squared") %>%
               unlist() %>%
               round(3),
             after = 7) %>%
  insert_row("Adj. R2",
             map(models, glance) %>%
               map(`[`, "adj.r.squared") %>%
               unlist() %>%
               round(3),
             after = 8) %>%
  insert_row("*** p < 0.001;  ** p < 0.01;  * p < 0.05. \nNote. Includes controls for the lagged non-proficiency rate in mathematics for all students and for white students, percent of teachers with less than 3 years of experience, percent of teachers with MAs or above, the percent of English Language Learners, the percent of students qualifying for free or reduced price lunch, and the natural log of total school enrollment. Heteroscedastic-robust standard errors clustered at the school level are in parentheses.",
             rep("", 6),
             after = 9) %>%
  insert_column(rep("", 10), after = 4) %>%
  set_colspan(10, 1, 8) %>%
  # 9 rows, 5 columns
  set_top_border(c(1, 3, 5), 1:8, 1) %>%
  set_top_border(2, c(2:4, 6:8), 1) %>%
  set_bottom_border(9, 1:8, 1) %>%
  set_all_padding(1:10, 1:8, 0) %>%
  set_align(1:10, 2:8, "center")

tab %>% quick_docx()

tab %>%
  as_flextable() %>%
  flextable::width(width = c(2.5, rep(1, 8)))