message("Beginning empirical replication...")
rm(list=ls())

library(bizicount)
suppressMessages(library(texreg))
suppressMessages(library(DHARMa))

# Load Data
dat = terror

# Scale covariates
dat[, c("xcoord", "ycoord", "pop","mtns")] = scale(dat[, c("xcoord", "ycoord", "pop", "mtns")])

# Marginal Poisson Formulas
fmla.ful = att.ful ~ pop + mtns + xcoord * ycoord
fmla.bok = att.bok ~ pop + mtns + xcoord * ycoord



# Bivariate Poisson
bivpois = bizicount(
  fmla.ful,
  fmla.bok,
  data = dat,
  cop = "frank",
  margins = c("pois", "pois"),
  keep = TRUE
)


# Use summary() to print standard errors with coefs
summary(bivpois)


# He, et al's test for zero modification on bivpois object
zi_test(bivpois)


# Define zero-inflated marginal formulas
fmla.zi.ful = att.ful ~ pop + mtns + xcoord * ycoord | 
  pop + mtns + xcoord * ycoord

fmla.zi.bok = att.bok ~ pop + mtns + xcoord * ycoord | 
  pop + mtns + xcoord * ycoord


# Bivariate, zero-inflated Poisson
zi_bivpois = bizicount(
  fmla.zi.ful,
  fmla.zi.bok,
  data = dat,
  cop = "frank",
  margins = c("zip", "zip"),
  keep = T
)


# Use summary() to print standard errors with coefs
summary(zi_bivpois)


# Make Dharma list for zero-inflated bizicount object (zi_bivpois)
zip_dharmas = make_DHARMa(zi_bivpois, nsim = 5000, seed = 789443)

class(zip_dharmas)
lapply(zip_dharmas, class)


# Test marginal residuals with DHARMa
lapply(zip_dharmas, testResiduals, plot = FALSE)



# Bivariate, zero-inflated negative-binomial
zi_bivnb = bizicount(
  fmla.zi.ful,
  fmla.zi.bok,
  data = dat,
  cop = "frank",
  margins = c("zinb", "zinb"),
  keep = T
)


# Dharma diagnostics on zero-inflated negative binomial
chi_stat_nb = -2*(logLik(zi_bivpois) - logLik(zi_bivnb))
pchisq(chi_stat_nb, df = 2, lower.tail = FALSE)

zinb_dharmas = make_DHARMa(zi_bivnb, nsim = 5000, seed = 12473)

lapply(zinb_dharmas, testResiduals, plot = FALSE)





# Texreg
mods = list(bivpois, zi_bivpois, zi_bivnb)

# Will produce latex code for table in manuscript,
texreg(mods,
       groups = list("Count Model" = 1:6, "Zero Inflation" = 7:12),
       reorder.coef = c(1:6, 8:13, 14, 7),
       digits = 2,
       custom.model.names = c(
         "Poisson: Fulani",
         "Poisson: BH",
         "ZIP: Fulani",
         "ZIP: BH",
         "ZINB: Fulani",
         "ZINB: BH"
       )
)

screenreg(mods,
         groups = list("Count Model" = 1:6, "Zero Inflation" = 7:12),
         reorder.coef = c(1:6, 8:13, 14, 7),
         digits = 2
)

# output to file for easy replication checking
screenreg(mods,
          groups = list("Count Model" = 1:6, "Zero Inflation" = 7:12),
          reorder.coef = c(1:6, 8:13, 14, 7),
          digits = 2,
          file= "Figures/table_3.txt"
)



## Appendix Dharma residual plots

name = c("fulani", "bh")

for(i in seq_along(name)){
  
  pdf(paste0("Figures/appendix_figure_B", i, ".pdf"), width=8, height=4.5)
    plot(zip_dharmas[[i]])
  dev.off()
    plot(zip_dharmas[[i]])
}


message("Empirical replication complete.")
