// Historical analysis (descriptive statistics and regression models)

use "cleaned_ALL.dta"

*********************
*********************
** desc statistics **
*********************
*********************

// for projection multipler computations
tab binary_e_vul year, m // 593
tab binary_e_par year, m // 1,106
tab binary_e_alg year, m // 777
tab binary_e_other year, m // 537
sum case_e_vul // 648
  return list
sum case_e_par // 1,856
  return list
sum case_e_alg // 914
  return list
sum case_e_other // 605
  return list

// geographic and temporal distribution of cases (exposure and reporting locations)
table state, contents(sum case_e_all sum case_r_all) column row
table year, contents(sum case_e_all sum case_r_all) column row
table state year, contents(sum case_e_all) column row
table state year, contents(sum case_r_all) column row
table state month, contents(sum case_e_all) column row
table state month, contents(sum case_r_all) column row

table year, contents(sum case_e_vul sum case_e_par sum case_e_alg sum case_e_other) column row
table region, contents(sum case_e_vul sum case_e_par sum case_e_other) column row
table region, contents(sum case_e_food sum case_e_non) column row

table month, contents(sum case_e_all mean sst) column row
table month if state!="HI" & state!="FL", contents(sum case_e_all) column row 

*lowess binary_e_all sst
*lowess binary_e_par sst
*lowess binary_e_vul sst
*lowess binary_e_other sst
*lowess binary_e_alg sst

*lowess binary_e_food sst
*lowess binary_e_non sst

*scatter sst sss
corr sst sss
*scatter sst tmean
corr sst tmean
*scatter tmean sss
corr sss tmean

*twoway scatter (sss sst), name(graph1, replace)
*twoway scatter (tmean sst), name(graph2, replace)
*graph combine graph1 graph2

*scatter case_e_all sst
*scatter case_e_all sst if case_e_all>0
*scatter case_e_all sss
*scatter case_e_all sss if case_e_all>0

*twoway scatter sst sss if binary_e_all==0, mcolor(navy) || scatter sst sss if binary_e_all==1, mcolor(cranberry)
  
table binary_e_all, contents(mean sst mean tmean mean sss)
table binary_e_vul, contents(mean sst mean tmean mean sss)
table binary_e_par, contents(mean sst mean tmean mean sss)

tabulate binary_e_all if month>=5 & month<=10

corr case_e_all case_r_all 
tabulate binary_e_all binary_r_all


*****************
*****************
** regressions **
*****************
*****************

************************************************************************************************
* logit, species, spline, PREFFERRED -- MAIN TEXT (TABLE 2), SUPPLEMENTAL MATERIALS (TABLE S9) *
************************************************************************************************

eststo clear 
// PARA 
mkspline sst1 8 sst2 = sst 
eststo: logit binary_e_par sst2 i.state_num ib(last).year, cluster(FIPS) 
drop sst1 sst2
// VUL
mkspline sst1 13 sst2 20 sst3 = sst 
eststo: logit binary_e_vul sst2 sst3 i.state_num ib(last).year, cluster(FIPS)
drop sst1 sst2 sst3
// ALG
mkspline sst1 5 sst2 22 sst3 = sst 
eststo: logit binary_e_alg sst2 sst3 i.state_num ib(last).year, cluster(FIPS)
drop sst1 sst2 sst3
// OTHER 
eststo: logit binary_e_other sst i.state_num ib(last).year, cluster(FIPS) 

*****************************************************************************
* logit, species, spline, with low end -- SUPPLEMENTAL MATERIALS (Table S5) *
*****************************************************************************

eststo clear 
// PARA 
mkspline sst1 8 sst2 = sst 
eststo: logit binary_e_par sst1 sst2 i.state_num ib(last).year, cluster(FIPS) 
drop sst1 sst2
// VUL
mkspline sst1 13 sst2 20 sst3 = sst 
eststo: logit binary_e_vul sst1 sst2 sst3 i.state_num ib(last).year, cluster(FIPS)
drop sst1 sst2 sst3
// ALG
mkspline sst1 5 sst2 22 sst3 = sst 
eststo: logit binary_e_alg sst1 sst2 sst3 i.state_num ib(last).year, cluster(FIPS)
drop sst1 sst2 sst3
// OTHER 
eststo: logit binary_e_other sst i.state_num ib(last).year, cluster(FIPS) 

****************************************************************************************************
* explore sensitivity of model to randomized exposure location - SUPPLEMENTAL MATERIALS (TABLE S6) *
******************************************************************* ********************************

// go back to "Cleaned_ExposureCounty.do", change seed, and re-run 
// seeds used:
//	98765 (USE THIS SEED THEREAFTER)
//	56789
//	12345
//	54321
//	11111
reg binary_e_all sst i.state_num i.year 

**********************************************************
* explore estimators - SUPPLEMENTAL MATERIALS (TABLE S7) *
**********************************************************

eststo clear 
eststo: reg binary_e_par sst i.state_num ib(last).year, cluster(FIPS)
eststo: logit binary_e_par sst i.state_num ib(last).year, cluster(FIPS)
margins, dydx(sst)
eststo: probit binary_e_par sst i.state_num ib(last).year, cluster(FIPS)
margins, dydx(sst)
eststo: reg binary_e_vul sst i.state_num ib(last).year, cluster(FIPS)
eststo: logit binary_e_vul sst i.state_num ib(last).year, cluster(FIPS)
margins, dydx(sst)
eststo: probit binary_e_vul sst i.state_num ib(last).year, cluster(FIPS)
margins, dydx(sst)

***********************************************************
* explore controls - SUPPLEMENTAL MATERIAL (Table S8) *
***********************************************************

eststo clear 
mkspline sst1 8 sst2 = sst 
eststo: logit binary_e_par sst2, cluster(FIPS) 
eststo: logit binary_e_par sst2 ib(last).year, cluster(FIPS) 
eststo: logit binary_e_par sst2 i.state_num, cluster(FIPS) 
eststo: logit binary_e_par sst2 i.state_num ib(last).year, cluster(FIPS) 
drop sst1 sst2

**********************************************************
* add salinity -- SUPPLEMENTAL MATERIALS (Table S10) *
**********************************************************

eststo clear 
// PARA 
mkspline sst1 8 sst2 = sst 
eststo: logit binary_e_par sst2 sss i.state_num ib(last).year if region=="gulf", cluster(FIPS) 
eststo: logit binary_e_par sst2 sss i.state_num ib(last).year if region=="east", cluster(FIPS) 
eststo: logit binary_e_par sst2 sss i.state_num ib(last).year if region=="west", cluster(FIPS) 
drop sst1 sst2
// VUL
mkspline sst1 13 sst2 20 sst3 = sst 
eststo: logit binary_e_vul sst2 sst3 sss i.state_num ib(last).year if region=="gulf", cluster(FIPS)
eststo: logit binary_e_vul sst2 sst3 sss i.state_num ib(last).year if region=="east", cluster(FIPS)
eststo: logit binary_e_vul sst2 sst3 sss i.state_num ib(last).year if region=="west", cluster(FIPS)
drop sst1 sst2 sst3
// ALG
mkspline sst1 5 sst2 22 sst3 = sst 
eststo: logit binary_e_alg sst2 sst3 sss i.state_num ib(last).year if region=="gulf", cluster(FIPS)
eststo: logit binary_e_alg sst2 sst3 sss i.state_num ib(last).year if region=="east", cluster(FIPS)
eststo: logit binary_e_alg sst2 sst3 sss i.state_num ib(last).year if region=="west", cluster(FIPS)
drop sst1 sst2 sst3
// OTHER
eststo: logit binary_e_other sst sss i.state_num ib(last).year if region=="gulf", cluster(FIPS) 
eststo: logit binary_e_other sst sss i.state_num ib(last).year if region=="east", cluster(FIPS) 
eststo: logit binary_e_other sst sss i.state_num ib(last).year if region=="west", cluster(FIPS)