***=======================================================*	 
* FEM scenario combination - lifetime
* Mar 2008
* May 2008
* Aug 2008, revised obesity scenarios
* Aug 30 2008, add cases avoided
***=======================================================*	 

clear
cap clear mata
set more off
set mem 500m

***=============================*	 
* DEFINE DIRECTORIES 
global workdir "/zeno/a/zyuhui/DOL/Model"
global indata "/zeno/a/zyuhui/DOL/Input"
global indata2 "/zeno/a/zyuhui/DOL/Indata"
global netindata "/homer/c/Retire/yuhui/rdata2"
global outdata "/zeno/a/zyuhui/DOL/Output"
global netoutdata "/homer/c/Retire/yuhui/simu"

* For local
global workdir "//zeno/zeno_a/zyuhui/ahg/Model_yh"
global indata "//zeno/zeno_a/zyuhui/ahg/Input_yh"
global indata2 "//zeno/zeno_a/zyuhui/ahg/Indata_yh"
global outdata "//zeno/zeno_a/zyuhui/ahg/Output_yh"

global netindata "//homer/homer_c/Retire/yuhui/rdata2"
global netoutdata "//homer/homer_c/Retire/yuhui/simu"


global workdir "//zeno/zeno_a/DOL/Othercodes_ahg"
global indata 	"//zeno/zeno_a/DOL/Input_yh"
global indata2 	"//zeno/zeno_a/DOL/Indata_yh"
global outdata 	"//zeno/zeno_a/DOL/Output_prev2"
global outdata 	"//zeno/zeno_a/DOL/Output_prev2"
global netindata "//homer/homer_c/Retire/yuhui/rdata2"
global netoutdata "//homer/homer_c/Retire/yuhui/simu"


#d;
global scnr status_quo Diabetes_R 
Diabetes_H Hbp_R Hbp_H Smoking_R Smoking_H Obese_R Obese_H Diabetes_E Hbp_E Smoking_E Obese_E;
#d cr

#d;
	global scnr status_quo 
	Diabetes_R Hbp_R Obs_R Smoking_R 
	Diabetes_J Hbp_J Obs_J Smoking_J
	Diabetes_E Hbp_E Obs_E Smoking_E	;
#d cr

#d; 
global grplist 5152_all 5152_hispan 5152_black 5152_white ;
global grplist 5152_all ;
#d cr

local disrate1 = 0
local disrate2 = 0.03
set trace off

* foreach fend in all hispan black white {
foreach fend in all {

	local f  5152
	local j = 0
	foreach x in $scnr {
		local j = `j' + 1
	}	
	
	#d;
	global clist SCRID LE HEALTH_LE DISABLE_LE QALY 
		GOVREV GOVEXP LFSPEND POP TREATED DIABE04 HIBPE04 DIABE_PREV HIBPE_PREV pop_lfdiabe pop_lfhibpe pop_lfsmoken pop_lfobese ; 
	#d cr
	
	local k = 0
	foreach x in $clist {
		local k = `k' + 1
	}

	#d;
	global sname status_quo none1 e25 DBT_25 HBP_25 OBS_25 SMK_25 
	 none2 e50 DBT_50 HBP_50 OBS_50 SMK_50
		none3 e100 DBT_100 HBP_100 OBS_100  SMK_100 ;
	#d cr
	
	#d;
	global sname status_quo 
		none1 Diabetes_R Hbp_R Obs_R Smoking_R 
		none2 Diabetes_J Hbp_J Obs_J Smoking_J
		none3 Diabetes_E Hbp_E Obs_E Smoking_E	;
	#d cr
	
	local j = 0
	foreach x in $sname {
		local j = `j' + 1
	}
		
	matrix matr = J(`j',`k',.)
	matrix rownames matr = $sname
	matrix colnames matr = $clist 
	matrix list matr
	
	local scrid = 0
	foreach scr in $sname {
		local scrid = `scrid' + 1
		matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"SCRID")] = `scrid'
			
		if !( strpos("`scr'", "none") > 0| "`scr'" == "e25" | "`scr'" == "e50" | "`scr'" == "e100" ) {
	
			drop _all
			use "$outdata//`f'_s`scr'_c.dta"
			
			gen govexp = ttl_ssben + ttl_ssiben + ttl_diben + ttl_caremd + ttl_caidmd
			gen govrev = ttl_ftax + ttl_stax + ttl_ctax + ttl_hoasi + ttl_hmed
			sort year, stable
			* Prevalence of diabetes at year 2004
			foreach d in diabe hibpe {
				gen `d'04 = `d'[1] 
			}
			* Total cases of diabetes and hiprevention
			foreach d in diabe hibpe {
				egen `d'_prev = total(`d'_new)
			}		
			* Total cases treated
			egen treated = total(treat_new)
			drop treat_new
			ren treated treat_new
			
			keep year pop qaly pop_healthy govrev govexp ttl_totmd treat_new diabe04 hibpe04 diabe_prev hibpe_prev pop_lfdiabe pop_lfhibpe pop_lfsmoken pop_lfobese
			
			* Include odds years
			qui sum year
			local myr = r(min)
			sort year, stable
			expand 2
			sort year, stable
			gen year_old = year
			replace year = `myr' + _n-1
			drop if year > 2054
			
			* For those already lifetime measure
			foreach x in treat_new diabe04 hibpe04 diabe_prev hibpe_prev {
				qui sum year
				replace `x' = 0 if year > r(min)
			}
			
			gen id = 1
			foreach x in pop qaly pop_healthy {
				gen `x'_old = `x'
				local k1 = _n+1
				local k2 = _n-1
				xtset id year
				sort id year
		
				replace `x' = (f.`x' + l.`x')/2 if _n > 1 & _n < _N & mod(year,2) == 1
				replace `x' = `x' * ((1/(1+`disrate1'))^(year-2004))
			}
		
			foreach x in govrev govexp ttl_totmd {
				gen `x'_old = `x'
				local k1 = _n+1
				local k2 = _n-1
				xtset id year
				sort id year
		
				replace `x' = (f.`x' + l.`x')/2 if _n > 1 & _n < _N & mod(year,2) == 1
				replace `x' = `x' * ((1/(1+`disrate2'))^(year-2004))
			}
	
			* LE, disable-free LE, lifetime ME
			
			qui sum pop if year == `myr'
			local init_pop = r(mean)
			qui sum pop
			local ttlyr = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"LE")] = `ttlyr'/`init_pop'
	
			qui sum qaly
			local ttlqaly = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"QALY")] = `ttlqaly'/`init_pop'
			
			qui sum pop_healthy
			local ttlyr = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"HEALTH_LE")] = `ttlyr'/`init_pop'	
			
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"DISABLE_LE")] = ///
			matr[rownumb(matr,"`scr'"), colnumb(matr,"LE")] - ///
			matr[rownumb(matr,"`scr'"), colnumb(matr,"HEALTH_LE")] 
			
			qui sum ttl_totmd 
			local ttlmd = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"LFSPEND")] = int(`ttlmd'*1000/`init_pop'	)
		
			qui sum govrev
			local ttl = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"GOVREV")] =int( `ttl'*1000/`init_pop')
		
			qui sum govexp
			local ttl = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"GOVEXP")] = int(`ttl'*1000/`init_pop')			
			
			qui sum pop if year == 2004
			local ttl_pop = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"POP")] = `ttl_pop'		
			
			dis "*************************"
			qui sum treat_new
			local ttl_treat = r(sum)
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"TREATED")] = `ttl_treat'/`ttl_pop'		
			
			qui sum year
			local minyear = r(min)
			qui sum diabe04 if year == `minyear'
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"DIABE04")] = r(mean)		
			qui sum hibpe04 if year == `minyear'
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"HIBPE04")] = r(mean)	
			qui sum diabe_prev
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"DIABE_PREV")] = r(sum)/`ttl_pop'
			qui sum hibpe_prev
			matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"HIBPE_PREV")] = r(sum)/`ttl_pop'
	
			/* Aug 30, 2008 */
			* For cases prevented, for every two years
			foreach x in pop_lfdiabe pop_lfhibpe pop_lfsmoken pop_lfobese {
			 	replace `x' = 0 if mod(year,2) != 0
				qui sum `x'
				matrix matr[rownumb(matr,"`scr'"), colnumb(matr,"`x'")] = r(sum)		
			}
		}
	}
	
	set trace off
	matrix matr_`fend' = matr
	matrix colnames matr_`fend' = $clist
	set trace off
}
drop _all

* foreach fend in all hispan black white {
foreach fend in all {
	svmat matr_`fend', names(col)
}

* Names of scenarios
gen scr = ""
local i = 1 
foreach s in $sname {
	sort SCRID, stable
	replace scr = "`s'" if SCRID == `i'
	local i = `i' + 1
}


	sort SCRID, stable
	foreach o in LE HEALTH_LE DISABLE_LE GOVREV GOVEXP LFSPEND POP QALY pop_lfdiabe pop_lfhibpe pop_lfobese pop_lfsmoken {
		qui sum `o'`fend' if SCRID == 1
		local base = r(mean)
		qui gen `o'`fend'_diff = (`o'`fend' - `base') if SCRID`grp' > 1
	}

	qui gen ECONGAIN_diff = QALY_diff * 100000 - LFSPEND_diff
	
* EFFECT OF TREATMENT ON THE TREATED
	foreach x in LE DISABLE_LE QALY LFSPEND ECONGAIN {
		gen `x'_treat = `x'_diff/TREATED
	}

* Proportion prevented
local totpop = 7.88
/*
	gen PREVENTED = pop_lfdiabe * 100 /10^6/`totpop' if strpos(scr, "DBT" ) > 0 & SCRID`grp' > 1
	replace PREVENTED = pop_lfhibpe *100 /10^6/`totpop'  if strpos(scr, "HBP" ) > 0 &  SCRID`grp' > 1	
	replace PREVENTED = pop_lfsmoken *100 /10^6/`totpop' if strpos(scr, "SMK") > 0 & SCRID`grp' > 1
	replace PREVENTED = pop_lfobese  *100 /10^6/`totpop'  if strpos(scr, "OBS") > 0 & SCRID`grp' > 1	 	
*/
	gen PREVENTED = pop_lfdiabe * 100 /10^6/`totpop' if strpos(scr, "Diabetes" ) > 0 & SCRID`grp' > 1
	replace PREVENTED = pop_lfhibpe *100 /10^6/`totpop'  if strpos(scr, "Hbp" ) > 0 &  SCRID`grp' > 1	
	replace PREVENTED = pop_lfsmoken *100 /10^6/`totpop' if strpos(scr, "Smoking") > 0 & SCRID`grp' > 1
	replace PREVENTED = pop_lfobese  *100 /10^6/`totpop'  if strpos(scr, "Obs") > 0 & SCRID`grp' > 1	 	

	#d;
	sort SCRID, stable;
	outsheet scr PREVENTED LE QALY  LFSPEND ECONGAIN_diff LE_treat DISABLE_LE_treat QALY_treat LFSPEND_treat ECONGAIN_treat using
	"$outdata/prevention_lifetime_v5.csv", replace comma nol ;
	#d cr