drop _all
capture log close
capture log using logit1.log, replace

*********************************************************
* Examples of logistic  regression models 
*********************************************************

use nhanesrh2
* look only at adults
drop if age < 18

//
more

***************************************
*       ANALYSIS PLAN
*          P   = dm
*  Key IndVars = age & female
*   Covariates = educ3 bmi sbp
***************************************
tab dm

//
more

global all "sbp age female educ3 bmi dm"

global x "age female i.educ3 bmi sbp"
global x2 "age i.educ3 bmi sbp"

//
more


* Model 1 (Sex)
logit dm female
//
more

di exp(-.0855)

logistic dm female
//
more

* Model 2 (Age)
logistic dm age
//
more

* Model 3 (Age & Sex)
logistic dm age female
//
more

* Full Model (Age & Sex plus all covariates)
xi: logistic dm $x

* Reformat 
//
more
reformat
//
more

*
* Discrimination (AUROC curve = C statistic)
* 
//
more
xi: logistic dm $x
lroc
graph export lroc1.emf, replace 

//
more
lroc, nog
//
more

*
* LR Test for improvements in model Discrimination
*   

qui xi: logit dm $x
keep if e(sample)
* So that we can compare models with the exact same observations
//
more
xi: logit dm age
lroc, nog
est store A
//
more
xi: logit dm age female
lroc, nog
lrtest A
//
more

xi: logit dm $x
lroc, nog
lrtest A
//
more

drop _all
use nhanesrh2
drop if age < 22

*
* Multicolinearity?
* 
//
more
*
* vif doesn't work with logit
* Must download collin.ado from web!!!
*

collin dm age female bmi sbp educ3

//
more

*
* Posterior Predictions using the predict & adjust commands
* 
//
more

xi: logistic dm age i.educ3
predict phat if e(sample)
drop if phat==.
sum dm phat
//
more

list age educ3 dm phat in 1/20

//
more

* What does the adjust command do?
//
more
xi: logistic dm age educ3
adjust age if e(sample), by(educ3) pr
//
more

adjust age = 40 if e(sample), by(educ3)  pr
adjust age = 60 if e(sample), by(educ3)  pr
adjust age = 80 if e(sample), by(educ3)  pr

//
more

* Review of Interaction (Conditional) Effects
//
more

* Stratified Analysis to see if age effects vary by bmi?
//
more

xi: logistic dm $x2 if bmi <30
xi: logistic dm $x2 if bmi >=30

//
more
* Interaction term to see if effect varies by sex?
recode bmi (10/29.9999999=0) (30/max=1), gen(bmi2)
gen agebmi = age*bmi2

//
more

logistic dm age bmi2 agebmi

//
more

adjust age = 30 agebmi=0 if e(sample) & bmi2==0, by(bmi2)  pr

adjust age = 30 agebmi=30 if e(sample) & bmi2==1, by(bmi2)  pr

adjust age = 70 agebmi=0 if e(sample) & bmi2==0, by(bmi2)  pr
adjust age = 70 agebmi=70 if e(sample) & bmi2==1, by(bmi2)   pr

//
more


*
* Hosmer-Lemeshow Tables
* 
//
more
sum dm phat

//
more

drop if phat == .
sort phat

li phat dm in 1/15

//
more
sort phat
gen cat5 = int(5*_n/(_N + 1))
bysort cat5: sum dm phat
*********************************************************************************
* Although 5 categories are used here, it is traditional to use deciles, IF
*  you have sufficient sample size to produce stable estimates for 10 subgroups
*********************************************************************************
//
more

graph bar (mean) phat dm, by(cat5, rows(1))


//
more

lfit, group(5)

//
more