/** Aggregate mortality across age and demographic groups, compute standardized mortality ratio and confidence interval
*/

log using logs/aggregate_hrsmort.log, replace

* confidence level for confidence intervals
local conflev 0.95

* Weir replication
tempfile weir
forvalues yr = 2001/2015 {
	append using data/weir_mortality_repl`yr'.dta
}
gen method = "Weir replication"
save `weir'

* Unweighted HRS death dates only
clear
tempfile hrsdeath
forvalues yr = 2001/2015 {
	append using data/hrsdeath_mortality`yr'.dta
}
gen method = "HRS (no NDI)"
save `hrsdeath'

* Weighted HRS death dates only
clear
tempfile weighted_hrsdeath
forvalues yr = 2001/2012 {
	append using data/weighted_hrsdeath_mortality`yr'.dta
}
gen method = "Weighted HRS (no NDI)"
save `weighted_hrsdeath'

* combine all methods into a single file
clear
append using `weir' `hrsdeath' `weighted_hrsdeath'
sort method year

* aggregate within each method and year
collapse (sum) exposure wtexposure died wtdied wt2died expdied wtexpdied, by(method year)

* compute mortality ratios
gen mortratio = died/expdied
gen wtmortratio = wtdied/wtexpdied


* compute confidence intervals

** compute upper and lower confidence interval bounds for unweighted mortality ratio using Ulm (1990; Chi-squared) method
gen X_L = invchi2(2*died,(1-`conflev')/2) / 2
gen X_U = invchi2(2*(died+1),1-(1-`conflev')/2) / 2
gen smr_lower = X_L / expdied
gen smr_upper = X_U / expdied

gen smr_lower_diff = mortratio - smr_lower
gen smr_upper_diff = smr_upper - mortratio

** compute upper and lower confidence interval bounds for weighted mortality ratio using Dobson et al (1991) method
gen T_L = wtdied + sqrt(wt2died / died) * (X_L - died)
gen T_U = wtdied + sqrt(wt2died / died) * (X_U - died)

gen wtsmr_lower = T_L / wtexpdied
gen wtsmr_upper = T_U / wtexpdied

gen wtsmr_lower_diff = wtmortratio - wtsmr_lower
gen wtsmr_upper_diff = wtsmr_upper - wtmortratio

drop X_L X_U T_L T_U

save data/mortality_ratios_allyrs.dta, replace

/*
restore

preserve

egen agegrp = cut(age), at(50(10)100)

collapse (sum) died expdied, by(year black male agegrp)

gen mortratio = died/expdied

drop if missing(black) | missing(male) | year > 2012

save hrs_mortality_allyrs_bydemog.dta, replace
*/
/*
restore

collapse (sum) wtexposed wtexposure wtdied, by(year black male age)
gen wtmortprobest = wtdied / wtexposed
gen wtmx = wtdied / wtexposure
gen wtqx = wtmx / (1 + wtmx/2)

sort year black male age

preserve
clear

forvalues yr = 2001/2015 {
	append using ../../base_data/cdc_nvsr_lifetable`yr'.dta
}

sort year black male age

tempfile lt
save `lt'

restore
merge 1:1 year black male age using `lt'
keep if _merge==3
drop _merge

li

keep year black male age qx wtqx
save qx_compare.dta, replace
*/


exit, clear