***********************************************************************
*                                                                     *
* framingham.do                                                       *
*                                                                     *
* PURPOSE:  analysis of the Framingham data with continuous           *
*           covariates                                                *
*           (1) specifying a functional relationship                  *
*           (2) evaluating fitted probabilities                       *
*           (3) summarizing accuracy                                  *
*                                                                     *
* DATE:  01/05/03                                                     *
*                                                                     *
***********************************************************************

***
*** input and label the data
***

#delimit ;
infile lexam 
       surv 
       cause 
       cexam 
       chd 
       cva 
       cancer 
       other 
       sex 
       age
       height
       weight
       chol1
       chol2
       dbp
       sbp
       mrw
       smoke    using fram-dat.txt;
#delimit cr

generate male = sex
recode male 2=0

generate bmi = (weight/2.2)/( (height/39)*(height/39) )

***
*** labels
***
label variable chd "death from CHD"
label variable cause "cause of death"
label variable height "height (inches)"
label variable weight "weight (pounds)"
label variable chol1 "serum cholesterol - exam 1"
label variable chol2 "serum cholesterol - exam 2"
label variable dbp "diastolic blood pressure"
label variable sbp "systolic blood pressure"
label variable smoke "cigarettes/day"

***
*** recode missing values
***

mvdecode height weight chol1 chol2 smoke, mv(-1)

*** 
*** subset for analysis = males 40+ w/o CHD at entry 
***

drop if sex>1 | age < 40 | cexam==1

*** drop subjects with missing values (not ideal, but...)

drop if chol1 > 1000 | smoke > 1000 | bmi > 1000 | bmi > 1000 

***
*** univariate summaries
***

summarize

table lexam
table surv
table cause
summarize age
summarize dbp 
summarize sbp
summarize smoke
summarize bmi
summarize chol1

***
*** generate categorical variables
***

generate Cage = age
recode Cage min/44.99=1 45/49.99=2 50/54.99=3 55/max=4
label define Lage 1 "< 45" 2 "45-50" 3 "50-55" 4 "55 +"
label values Cage Lage

generate Cdbp = dbp
recode Cdbp min/79.99=1 79/99.99=2 100/119.99=3 120/max=4
label define Ldbp 1 "50-80" 2 "80-100" 3 "100-120" 4 "120+"
label values Cdbp Ldbp

generate Csbp = sbp
recode Csbp min/119.99=1 120/139.99=2 140/159.99=3 160/199.99=4 200/max=5
label define Lsbp 1 "<120" 2 "120-140" 3 "140-160" 4 "160-200" 5 "200+"
label values Csbp Lsbp

generate Csmoke = smoke
recode Csmoke 0=1 1/15=2 16/30=3 30/max=4
label define Lsmoke 1 "none" 2 "1-15" 3 "16-30" 4 "30+"
label values Csmoke Lsmoke

generate Cbmi = bmi
recode Cbmi min/21.999=1 22/25.999=2 26/29.999=3 30/max=4
label define Lbmi 1 "<22" 2 "22-26" 3 "26-30" 4 "30+"
label values Cbmi Lbmi

generate Cchol1 = chol1
recode Cchol1 min/200=1 201/250=2 251/300=3 301/max=4
label define Lchol1 1 "<=200" 2 "201-250" 3 "251-300" 4 "300+"
label values Cchol1 Lchol1

***
*** univariate analysis with CHD
***

*--------------------------------------------------------------------------
*******
* (1) *  BMI
*******

table Cbmi
table Cbmi chd
tabodds chd Cbmi, or

xi: logistic chd i.Cbmi
logit

predict pfit1
predict logit1, xb

graph logit1 bmi
graph pfit1 bmi

********************
*** linear bmi   ***
********************
generate bmi2c = (bmi-25)*(bmi-25)
generate bmi1c = (bmi-25)

logistic chd bmi1c
logit
lrtest, saving(2)

predict pfit2
predict logit2, xb
graph logit2 bmi
graph pfit2 bmi

***********************
*** quadratic bmi   ***
***********************
logistic chd bmi1c bmi2c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3
predict logit3, xb
graph logit3 bmi
graph pfit3 bmi

***********************
*** square root bmi ***
***********************
generate sqrtbmi = sqrt(bmi)

logistic chd bmi sqrtbmi
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4
predict logit4, xb
graph logit4 bmi
graph pfit4 bmi

***************************
*** linear spline bmi   ***  (optional!)
***************************

mkspline bmiS1 24 bmiS2 28 bmiS3 = bmi

logistic chd bmiS1-bmiS3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5
predict logit5, xb
graph logit5 bmi, xline(24 28)
graph pfit5 bmi, xline(24 28)

***********************************
*** fractional polynomial bmi   ***  (optional!)
***********************************

fracpoly logistic chd bmi
logit

predict pfit6
predict logit6, xb
graph logit6 bmi
graph pfit6 bmi

fracplot bmi

*--------------------------------------------------------------------------
*******
* (2) *  AGE
*******

table Cage
table Cage chd
tabodds chd Cage, or

xi: logistic chd i.Cage
logit

predict pfit1b
predict logit1b, xb

graph logit1b age

********************
*** linear age   ***
********************
generate age2c = (age-50)*(age-50)
generate age1c = (age-50)

logistic chd age1c
logit
lrtest, saving(2)

predict pfit2b
predict logit2b, xb
graph logit2b age

***********************
*** quadratic age   ***
***********************
logistic chd age1c age2c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3b
predict logit3b, xb
graph logit3b age

***********************
*** square root age ***
***********************
generate sqrtage = sqrt(age)

logistic chd age sqrtage
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4b
predict logit4b, xb
graph logit4b age

***************************
*** linear spline age   ***  (optional!)
***************************

mkspline ageS1 45 ageS2 55 ageS3 = age

logistic chd ageS1-ageS3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5b
predict logit5b, xb
graph logit5b age, xline(45 55)

***********************************
*** fractional polynomial age   ***  (optional!)
***********************************

fracpoly logistic chd age
logit

predict pfit6b
predict logit6b, xb
graph logit6b age

fracplot age

*--------------------------------------------------------------------------
*******
* (3) *  SMOKE
*******

table Csmoke
table Csmoke chd
tabodds chd Csmoke, or

xi: logistic chd i.Csmoke
logit

predict pfit1c
predict logit1c, xb

graph logit1c smoke

********************
*** linear smoke   ***
********************
generate smoke2c = (smoke-10)*(smoke-10)
generate smoke1c = (smoke-10)

logistic chd smoke1c
logit
lrtest, saving(2)

predict pfit2c
predict logit2c, xb
graph logit2c smoke

***********************
*** quadratic smoke   ***
***********************
logistic chd smoke1c smoke2c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3c
predict logit3c, xb
graph logit3c smoke

***********************
*** square root smoke ***
***********************
generate sqrtsmk = sqrt(smoke)

logistic chd smoke sqrtsmk
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4c
predict logit4c, xb
graph logit4c smoke

***************************
*** linear spline smoke   ***  (optional!)
***************************

mkspline smokeS1 10 smokeS2 25 smokeS3 = smoke

logistic chd smokeS1-smokeS3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5c
predict logit5c, xb
graph logit5c smoke, xline(10 25)

***********************************
*** fractional polynomial smoke   ***  (optional!)
***********************************

fracpoly logistic chd smoke
logit

predict pfit6c
predict logit6c, xb
graph logit6c smoke

fracplot smoke

*--------------------------------------------------------------------------
*******
* (4) *  CHOL1
*******

table Cchol1
table Cchol1 chd
tabodds chd Cchol1, or

xi: logistic chd i.Cchol1
logit

predict pfit1d
predict logit1d, xb

graph logit1d chol1

********************
*** linear chol1   ***
********************
generate chol12c = (chol1-200)*(chol1-200)
generate chol11c = (chol1-200)

logistic chd chol11c
logit
lrtest, saving(2)

predict pfit2d
predict logit2d, xb
graph logit2d chol1

***********************
*** quadratic chol1   ***
***********************
logistic chd chol11c chol12c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3d
predict logit3d, xb
graph logit3d chol1

***********************
*** square root chol1 ***
***********************
generate sqrtchol = sqrt(chol1)

logistic chd chol1 sqrtchol
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4d
predict logit4d, xb
graph logit4d chol1

***************************
*** linear spline chol1   ***  (optional!)
***************************

mkspline chol1S1 200 chol1S2 275 chol1S3 = chol1

logistic chd chol1S1-chol1S3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5d
predict logit5d, xb
graph logit5d chol1, xline(200 275)

***********************************
*** fractional polynomial chol1   ***  (optional!)
***********************************

fracpoly logistic chd chol1
logit

predict pfit6d
predict logit6d, xb
graph logit6d chol1

fracplot chol1


*--------------------------------------------------------------------------
*******
* (5) *  DBP
*******

table Cdbp
table Cdbp chd
tabodds chd Cdbp, or

xi: logistic chd i.Cdbp
logit

predict pfit1e
predict logit1e, xb

graph logit1e dbp

********************
*** linear dbp   ***
********************
generate dbp2c = (dbp-100)*(dbp-100)
generate dbp1c = (dbp-100)

logistic chd dbp1c
logit
lrtest, saving(2)

predict pfit2e
predict logit2e, xb
graph logit2e dbp

***********************
*** quadratic dbp   ***
***********************
logistic chd dbp1c dbp2c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3e
predict logit3e, xb
graph logit3e dbp

***********************
*** square root dbp ***
***********************
generate sqrtdbp = sqrt(dbp)

logistic chd dbp sqrtdbp
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4e
predict logit4e, xb
graph logit4e dbp

***************************
*** linear spline dbp   ***  (optional!)
***************************

mkspline dbpS1 80 dbpS2 100 dbpS3 = dbp

logistic chd dbpS1-dbpS3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5e
predict logit5e, xb
graph logit5e dbp, xline(80 100)

***********************************
*** fractional polynomial dbp   ***  (optional!)
***********************************

fracpoly logistic chd dbp
logit

predict pfit6e
predict logit6e, xb
graph logit6e dbp

fracplot dbp

*--------------------------------------------------------------------------
*******
* (6) *  SBP
*******

table Csbp
table Csbp chd
tabodds chd Csbp, or

xi: logistic chd i.Csbp
logit

predict pfit1f
predict logit1f, xb

graph logit1f sbp

********************
*** linear sbp   ***
********************
generate sbp2c = (sbp-120)*(sbp-120)
generate sbp1c = (sbp-120)

logistic chd sbp1c
logit
lrtest, saving(2)

predict pfit2f
predict logit2f, xb
graph logit2f sbp

***********************
*** quadratic sbp   ***
***********************
logistic chd sbp1c sbp2c
logit
lrtest, saving(3)

lrtest, using(3) model(2)

predict pfit3f
predict logit3f, xb
graph logit3f sbp

***********************
*** square root sbp ***
***********************
generate sqrtsbp = sqrt(sbp)

logistic chd sbp sqrtsbp
logit
lrtest, saving(4)

lrtest, using(4) model(2)

predict pfit4f
predict logit4f, xb
graph logit4f sbp

***************************
*** linear spline sbp   ***  (optional!)
***************************

mkspline sbpS1 120 sbpS2 150 sbpS3 = sbp

logistic chd sbpS1-sbpS3
logit
lrtest, saving(5)

lrtest, using(5) model(2)

predict pfit5f
predict logit5f, xb
graph logit5f sbp, xline(120 150)

***********************************
*** fractional polynomial sbp   ***  (optional!)
***********************************

fracpoly logistic chd sbp
logit

predict pfit6f
predict logit6f, xb
graph logit6f sbp

fracplot sbp

*--------------------------------------------------------------------------