drop _all
capture log close
set more 1
capture log using NHANES_logit.log, replace

*********************************************************
*      Data Analysis  - Sample Stata Program     
*			NHANES_logit.do                       *
*    Evaluation of Predictors for Health Status 	  *
* Bivariate & Multivariate Results for DepVar=healthst  *
*        For CSP Summer Seminar in 2004                 *
*********************************************************

*Things to keep in mind regarding logistic regression:
*--sample size
*--chi2 test for model
*--what the coefficients mean

*"Goodness of Fit"
*--discrimination/predictiveness (AUROC)
*--callibration

use NHANES_transformed
* For DepVars = healthst
* Key IndVars = age, gender, race/ethnicity & educ

* look only at adults
drop if age<18

quietly {

* Make a global macro command for full casemix model
global demog "age female i.race"
global social "i.educ3 i.income"
global health "bmi"
global all "$demog $social $health"
}
*
*  For Table 3:   Results for DepVar = poorhlth(dichotomous DepVar)
*
gen poorhlth=healthst 
recode poorhlth 1/3=0 4/5=1 
label define poorhlth 0 "good+ health" 1 "fair/poor health"
label values poorhlth poorhlth 

tab poorhlth 


tab healthst poorhlth 
tab poorhlth, sum(age)

logistic poorhlth age

**what if we look at age by decade? (highest age cat is 85+)
gen age_dec=age
recode age_dec 18/19=1 20/29=2 30/39=3 40/49=4 50/59=5 60/69=6 70/79=7 80/85=8
tab age_dec,sum(age)

logistic poorhlth age_dec

**generate a dichotomous variable for race=black
gen AA=race
recode AA 1 3 4 =0 2=1

for var AA educ3 : tab poorhlth X, co ro ch

**How much statistical power do you have to have to detect differences in healthst by race?

sampsi .23 .28,p(.8)
sampsi .23 .28,p(.8) r(.065)
sampsi .23 .30,p(.8) r(.065)

**Get adjusted and unadjusted Odds Ratios
logistic poorhlth age 
lroc
graph export lroc1.tif, replace

logistic poorhlth female 
lroc
graph export lroc2.tif, replace 

qui xi:logistic poorhlth i.race 
reformat
lroc
graph export lroc3.tif, replace 

**Controlling for demographics
qui xi: logistic poorhlth $demog
reformat
lroc
graph export lroc4.tif, replace 

**and adding sociodemographics
qui xi: logistic poorhlth $demog $social
reformat
lroc
graph export lroc5.tif, replace 

**controlling for full casemix model
qui xi: logistic poorhlth $all
reformat
lroc
graph export lroc6.tif, replace 

* Making Predicted Probabilities 
predict phat

sum poorhlth phat
sum poorhlth if phat~=.
sum poorhlth if phat==.

* Looking at & testing for Calibration and Predcitiveness, in this example we're * creating 6 categories and here we're doing it all by hand
drop if phat==.  |  poorhlth ==.
sort phat
gen cat6 = int(6*_n/(_N + 1))
sort cat6
by cat6:sum poorhlth phat
graph bar (mean) poorhlth phat, by(cat6, rows(1))
graph export hl_bar1.tif, replace 

* Or we can just use this Hosmer-Lemshow command
hl2, group(6) plot2