**************************
* Match MEPS functional limitation measures with that in HRS
* Deal with IADL separately
* Calculate eq-5D and SF-6D in MEPS
**************************
clear

set more off
set mem 500m

global wkdir "."
global mepsdir "."
global hrsdir "."

cap log close
log using "$wkdir\eq5d_MEPS_HRS.log", replace

	global cov1 "male black hispan hsdrop somecol collgrad regnth regmid regwst widowed single"
	global cov2 "cancre diabe hibpe hearte lunge stroke depress"

	***********************************************
	* MEPS
	/*
	311 Depressive disorder, not elsewhere classified
	Depressive disorder NOS
	Depressive state NOS
	Depression NOS
	Excludes:
	acute reaction to major stress with depressive symptoms (308.0)
	affective personality disorder (301.10-301.13)
	affective psychoses (296.0-296.9)
	brief depressive reaction (309.0)
	depressive states associated with stressful events (309.0-309.1)
	disturbance of emotions specific to childhood and adolescence, with misery and unhappiness (313.1)
	mixed adjustment reaction with depressive symptoms (309.4)
	neurotic depression (300.4)
	prolonged depressive adjustment reaction (309.1)
	psychogenic depressive psychosis (298.0)
	*/
	***********************************************
	
	global deplist 
	
	***********************************************
	*Depression condition from condition file
	***********************************************
	use "$mepsdir/condition2001.dta", clear
	gen depress = 0
	foreach c in 296 311 {
		replace depress = 1 if inlist(icd9codx,"`c'")
	}
	tab icd9codx if depress == 1
	keep if depress == 1
	bys dupersid: keep if _n == 1
	keep dupersid depress
	sort dupersid
	save tmp,replace
	
	use "$mepsdir\MEPS9904_selected.dta"
	keep if yr == 2001
	sort dupersid
	merge dupersid using tmp
	tab _merge
	drop if _merge == 2
	erase tmp.dta
	replace depress = 0 if _merge != 3
	drop _merge

	
	* Replace the lung disease variable using ICD codes  
	replace lunge = lungecr 
	ren cancrecr cancre 
	gen age6574 = inrange(age,65,74)
	gen age75p  = age >= 75
	gen age70p  = age >= 70
	sort dupersid, stable
	save tmp, replace

	drop _all
	use "$mepsdir\csd2001.dta"
	merge dupersid using tmp, sort
	erase tmp.dta
	tab _merge
	keep if _merge == 3
	drop _merge
		
	* Keep aged 51 and over
	keep if age >= 51

	********************************************
	* Construct EQ-5D according to Shaw JW et al. 2005
	********************************************	
	gen mo = admobi42 if admobi42 > 0
	gen sc = adself42 if adself42 > 0
	gen ua = adacti42 if adacti42 > 0
	gen pd = adpayn42 if adpayn42 > 0
	gen ad = addepr42 if addepr42 > 0
	foreach x in mo sc ua pd ad {
		ren `x' `=upper("`x'")'
	}
	  
	/* Generate Dummy Variables for Levels 2 and 3 of Five Dimensions */
	gen m1 =0
	gen m2 =0
	gen s1 =0
	gen s2 =0
	gen u1 =0
	gen u2 =0
	gen p1 =0
	gen p2 =0
	gen a1 =0
	gen a2 =0
	
	replace m1 = 1 if MO == 2
	replace m2 = 1 if MO == 3
	replace s1 = 1 if SC == 2
	replace s2 = 1 if SC == 3
	replace u1 = 1 if UA == 2
	replace u2 = 1 if UA == 3
	replace p1 = 1 if PD == 2
	replace p2 = 1 if PD == 3
	replace a1 = 1 if AD == 2
	replace a2 = 1 if AD == 3
	
	
	/* Generate Interaction Terms (I2, I2-squared, I3, I3-squared) */
	gen m0 = 1 if m1 == 0 & m2 == 0
	gen s0 = 1 if s1 == 0 & s2 == 0
	gen u0 = 1 if u1 == 0 & u2 == 0
	gen p0 = 1 if p1 == 0 & p2 == 0
	gen a0 = 1 if a1 == 0 & a2 == 0
	
	egen i2 = neqany(m1 s1 u1 p1 a1), values(1)
	replace i2 = i2 - 1
	replace i2 = 0 if i2<0
	gen i22 = i2^2
	
	egen i3 = neqany(m2 s2 u2 p2 a2), values(1)
	replace i3 = i3 - 1
	replace i3 = 0 if i3<0
	gen i32 = i3^2
	
	
	/* Generate D1 Term */
	egen i1 = neqany(m0 s0 u0 p0 a0), values(1)
	gen d1 = 4 - i1
	replace d1 = 0 if d1<0
	
	

	/* Generate Raw Index Score */
	gen pred = .146016*m1 + .557685*m2 + .1753425*s1 + .4711896*s2 + .1397295*u1 + .3742594*u2 + .1728907*p1 + .5371011*p2 + .156223*a1 + .4501876*a2 + -.1395949*d1 + .0106868*i22  + -.1215579*i3 + -.0147963*i32
	gen eq5d = 1 - pred
	replace eq5d = . if MO == . | SC == . | UA == . | PD == . | AD == . 
	
	
	/* Drop Variables Generated by Program */
	drop m1-pred 
	
	/* Examine distribution */
	label var eq5d "EQ-5D score with US TTO weighting, Shaw JW Med Care 2005,43(3) 203-220"
	
	********************************************
	* Recode functional limitation variables
	********************************************
	cap drop iadl1 adl1p
	gen iadl1 = 0
	gen adl1p = 0
	foreach r in 53 {
		replace iadl1 = 1 if iadlhp`r' == 1
		replace adl1p = 1 if adlhlp`r' == 1
	}
	foreach r in 53 {
		replace iadl1 = . if iadl1 == 0 & iadlhp`r' < 0 
		replace adl1p = . if adl1p == 0 & adlhlp`r' < 0 
	}
	
	* Major activity limitation - work/housework
	gen mjlim = (actlim53 == 1)|(actlim53 == 1)  if inlist(actlim53,1,2)
	
	* Unable to carry out major activities
	gen mjuab = unable53 == 1 if mjlim == 1
	replace mjuab = 0 if mjlim == 0
		
	* Other any limitation
	gen anylim = wlklim53 == 1 | seedif42 == 1 | heardi42 == 1
	replace anylim = . if anylim == 0 & (!inlist(wlklim53,1,2) | !inlist(seedif42,1,2) | !inlist(heardi42,1,2))
	
	* Funcstat
	gen funclim = mjlim == 1 | anylim == 1 | adl1p == 1
	replace funclim = . if funclim == 0 & (mjlim == . | anylim == . | adl1p == .)
	gen iadlim = iadl1
	gen iadlim2 = iadlim
	replace iadlim2 = 0 if funclim == 1

	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = perwt]
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = perwt], stats(sd)

	replace iadl1 = 0 if adl1p == 1
	replace mjuab = 0 if iadl1 == 1 | adl1p == 1
	replace mjlim = 0 if mjuab == 1 | iadl1 == 1 | adl1p == 1
	replace anylim = 0 if mjlim == 1 | mjuab == 1 | iadl1 == 1 | adl1p == 1
	
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = perwt]
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = perwt], stats(sd)
	
	table funclim iadlim [aw = perwt], c(mean eq5d sd eq5d freq)
	tab funclim iadlim [aw = perwt], cell nofreq
	saveold "$mepsdir\MEPS_EQ5D.dta", replace	
	
	********************************************
	* Relationship between EQ-5D and functional status and disease conditions
	********************************************
	
	gen func_iadl = funclim * iadlim
	label var func_iadl "Activity limitation * IADL limitation"
	label var funclim "Activity limitation(ADL/Vision/Hearing/other)"
	label var iadlim "Any IADL limitation"
	
	foreach v in  eq5d  {
		reg  `v' funclim iadlim func_iadl [aw = perwt]
		est store e`v'_1
		reg  `v' iadl1 adl1p [aw = perwt]	
		est store e`v'_2
		reg `v'  $cov2  [pw = perwt]
		est store e`v'_3
		reg `v'  funclim iadlim func_iadl $cov2 [pw = perwt]
		est store e`v'_4
		matrix m`v' = e(b)
		reg `v'  iadl1 adl1p $cov2 [pw = perwt]
		est store e`v'_5			
		matrix m`v'_old = e(b)
		reg  `v' age6574 age75p $cov1 $cov2 funclim iadlim func_iadl [aw = perwt]
		est store e`v'_6
	}

	#d;
	estout eeq5d_* 
 	using "$wkdir\EQ5D2001.txt", 
		  cells(b(star fmt(%9.3f)) se(par fmt(%9.3f)))
              stats(r2_a N, fmt(%9.2g %9.0fc) label("Adjusted R2" "N"))
		  starlevels(* 0.10 ** 0.05 *** 0.01)
		  legend label collabel(,none)
              varlabels(_cons "Constant")
              prefoot("")
		  postfoot("")
              varwidth(10)
		  modelwidth(20)
		  replace;
	#d cr
	global cov3 "funclim iadlim mjlim anylim iadl1 adl1p"
	collapse $outvar  age6574 age75p $cov2 $cov3 [aw = perwt]
	gen data = "MEPS 2001"
	save tmp_meps.dta, replace
	drop _all

	*******************
	*HRS
	*******************	
	*Depression based on CESD score
	use hhidpn r*cesd using "$hrsdir/rndhrs_i.dta", clear
	forvalues i = 2/8 {
		ren r`i'cesd cesd`i'
	}
	reshape long cesd, i(hhidpn) j(wave)
	keep if wave == 4
	replace cesd = cesd >= 3 if cesd < . 
	ren cesd depress
	sort hhidpn wave
	save hrs_depress, replace
	

	* Year 1998
	drop _all
	use "$hrsdir\hrs_long_processed.dta"
	keep if age >= 51 & wave == 4 & iwstat == 1
	sort hhidpn wave
	merge hhidpn wave using hrs_depress
	tab _merge
	keep if _merge == 3
	drop _merge
	sort hhidpn 
	erase hrs_depress.dta
	save tmp,replace

	drop _all
	use "$hrsdir\h98f2b.dta"
	keep hhidpn f4017 f4021 f4022 f4032 f4201 f4205 f4206 f4238
	merge hhidpn using tmp,sort
	tab _merge
	keep if _merge == 3
	drop _merge
	erase tmp.dta
	* Combine working, housework, and any other limitations
	/*
	f4017           byte   %8.0g                  gd1. health prob
	f4021           byte   %8.0g                  gd1d.limit housework
	f4022           byte   %8.0g                  gd1e.limit in anyway
	f4032           byte   %8.0g                  gd2. keep from working
	f4201           byte   %8.0g                  gj1.health prob
	f4205           byte   %8.0g                  gj1d.limit housework
	f4206           byte   %8.0g                  gj1e.limit in any way
	f4238           byte   %8.0g                  gj13.keep from working
	*/
	
	cap drop iadl1
	cap drop adl1p
	gen iadl1 = iadla >= 1 if !missing(iadla)
	gen adl1p = adla >= 1 if !missing(adla)
	
	gen mjlim = inlist(f4017,1) | inlist(f4021,1) | inlist(f4201,1 ) | inlist(f4205,1)
	replace mjlim = . if mjlim == 0 & ((missing(f4017) & missing(f4201)) | (missing(f4021) & missing(f4205)))
	gen anylim = mjlim == 1 | f4022 == 1 | f4206 == 1
	replace anylim = . if (missing(mjlim) | (missing(f4022)& missing(f4206))) & anylim == 0
	
	* Don't make IADL1p and ADL lump together
	gen funclim = adl1p == 1 | mjlim == 1 | anylim == 1
	replace funclim =. if funclim == 0 & (adl1p == . | mjlim == . | anylim == .)
	
	keep if age >= 51
	gen iadlim = iadl1
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = wtresp]
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = wtresp], stats(sd)
	
	replace iadl1 = 0 if adl1p == 1
	replace mjlim = 0 if iadl1 == 1 | adl1p == 1
	replace anylim = 0 if mjlim == 1 | iadl1 == 1 | adl1p == 1
	
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = wtresp]
	tabstat funclim iadlim mjlim anylim iadl1 adl1p [aw = wtresp], stats(sd)

	* Nursing home
	tab nrshom funclim
	replace funclim = 1 if nrshom == 1

	replace adl1p = 1 if nrshom == 1
	replace iadl1 = 0 if nrshom == 1
	gen func_iadl = funclim * iadlim
	
	* local a = meq_index[1, colnumb(meq_index,"adl1p")]
	* matrix ad = (`a'-0.1)
	* matrix ad = (-0.1)
	* matrix colnames ad = "nrshom"
	
	* matrix meq5d = meq5d, ad
	matrix score eq5d_p = meq5d
	* matrix score eq5d_pold = meq5d_old
	* matrix score eq5d_pL = meq5dL
	* matrix score eq5d_poldL = meq5dL_old

	table funclim iadlim, c(mean eq5d_p sd eq5d_p freq)
	tab funclim iadlim [aw = wtresp], cell nofreq
	
	reg eq5d_p cancre diabe hearte hibpe lunge stroke iadl1 adl12 adl3 depress [aw = wtresp]
	est store eq5d_hrs
	
	* Store the data for ZENO estimation
	keep eq5d_p cancre diabe hearte hibpe lunge stroke iadl1 adl12 adl3  wtresp depress
	save "$hrsdir\EQ5D_HRS.dta", replace
		
	#d;
	estout eq5d_hrs
 	using "$wkdir\EQ5D_hrs.txt", 
		  cells(b(star fmt(%9.3f)) se(par fmt(%9.3f)))
              stats(r2_a N, fmt(%9.2g %9.0fc) label("Adjusted R2" "N"))
		  starlevels(* 0.10 ** 0.05 *** 0.01)
		  legend label collabel(,none)
              varlabels(_cons "Constant")
              prefoot("")
		  postfoot("")
              varwidth(10)
		  modelwidth(20)
		  replace;
	#d cr
	

log close