**** Analyses replicating results for Rehill et at 2013 Submitted. ***
insheet using london_annual_analysis_data.csv, comma case clear

* create putative data break indicators and year spline
foreach year in 1957 1958 1964 1965 1966 1967 1968 1983 1984 1974 1975 {
 gen indic`year'=(ws_year>=`year') 
 }
mkspline ncs6year=ws_year, cubic nk(6) 

**************************************** Table 1 ***********************************************
glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
glm AGE65PLUS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
glm AGEu65    indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
glm CVD       indic1965 indic1966  indic1967 indic1968 indic1983 indic1984 ///
                                   propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
glm RESP      indic1965 indic1966  indic1967 indic1968 indic1983 indic1984 ///
                                   propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef

 
 **************************** SENSITIVITY ANALYSES ****************************
* Baseline: main model
 glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
 
* Alternative count distribution assumptions 
 glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, link(log)  // slr
 nbreg ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat  // slr

* Alternative time splines
 mkspline ncs3year=ws_year, cubic nk(3) 
 glm ALLDEATHS indic1965 indic1966  propinf ncs3year* annual_cold annual_heat, f(p) sca(x2) ef
 
 spbase ws_year, nk(7) gen( ncs9year) // mkspline only runs up to 7 total knots = 6df
 gen ncs9year0=ws_year
 glm ALLDEATHS indic1965 indic1966  propinf ncs9year* annual_cold annual_heat, f(p) sca(x2) ef

*Additional steps
 glm ALLDEATHS indic1965 indic1966 indic1957 indic1958 propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
 glm ALLDEATHS indic1965 indic1966 indic1964  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef
 glm ALLDEATHS indic1965 indic1966 indic1974 indic1975 propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef

* Alternative cold and heat thresholds
 glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold_t15 annual_heat_t15, f(p) sca(x2) ef
 glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold_t21 annual_heat_t21, f(p) sca(x2) ef
 glm ALLDEATHS indic1965 indic1966  propinf ncs6year* annual_cold annual_heat, f(p) sca(x2) ef


* No adj for flu
 glm ALLDEATHS indic1965 indic1966 ncs6year* annual_cold annual_heat, f(p) sca(x2) ef

***************************************************************************

 
********************** RESIDUAL ANALYSIS **********************************
 gen logall=log(ALLDEATHS)
 glm logall indic1965 indic1966  propinf ncs6year* annual_cold annual_heat  
 scalar sdres=sqrt(e(dispers))
 predict resid, res
 predict fv
 gen stresid=resid/sdres
 summ  stresid, detail
 twoway scatter stresid ws_year
 twoway scatter stresid fv, xtitle("Fitted log count")
 sort ws_year
 tsset ws_year 
 pac stres
 ***************************************************************************