cd "/Users/steve/Documents/teaching/c2016/biostat208/labs/lab3/"
* do "/Users/steve/Documents/teaching/c2016/biostat208/labs/lab3/lab3.do"

log using lab3, replace
use lab3, clear

tab physact, sum(bmi)
reg bmi i.physact
* again, showing reference level
reg bmi i.physact, baselevels
* compute mean BMI at each level of physact
forvalues i = 1/5 {
lincom _cons + `i'.physact
}
* again, using margins
margins physact

* test for heterogeneity
testparm i.physact
* test for linear trend
contrast q(1).physact, noeffects
* test for departure from linear trend
contrast q(2/4).physact, noeffects

* repeat with adjustment for age, smoking, and alcohol use
reg bmi i.physact age10 smoking drnkspwk
* compute marginal mean BMI at each level of physact, 
* averaging over other predictors
margins physact
* verifiy results from last command
gen y1 = _b[_cons] + _b[1.physact]*1 + _b[age10 ]*age10 + _b[smoking]*smoking + _b[drnkspwk]*drnkspwk if e(sample)
gen y2 = _b[_cons] + _b[2.physact]*1 + _b[age10 ]*age10 + _b[smoking]*smoking + _b[drnkspwk]*drnkspwk if e(sample)
gen y3 = _b[_cons] + _b[3.physact]*1 + _b[age10 ]*age10 + _b[smoking]*smoking + _b[drnkspwk]*drnkspwk if e(sample) 
gen y4 = _b[_cons] + _b[4.physact]*1 + _b[age10 ]*age10 + _b[smoking]*smoking + _b[drnkspwk]*drnkspwk if e(sample) 
gen y5 = _b[_cons] + _b[5.physact]*1 + _b[age10 ]*age10 + _b[smoking]*smoking + _b[drnkspwk]*drnkspwk if e(sample) 
summ y1 y2 y3 y4 y5

* compute marginal mean BMI at each level of physact, 
* with other predictors fixed at mean values
margins physact, atmeans
* verifiy results from last command
* first generate mean values for predictors in estimation sample
egen m_age10 = mean(age10) if e(sample)
egen m_smoking = mean(smoking) if e(sample)
egen m_drnkspwk = mean(drnkspwk) if e(sample)
gen z1 = _b[_cons] + _b[1.physact]*1 + _b[age10 ]*m_age10 + _b[smoking]*m_smoking + _b[drnkspwk]*m_drnkspwk if e(sample)
gen z2 = _b[_cons] + _b[2.physact]*1 + _b[age10 ]*m_age10 + _b[smoking]*m_smoking + _b[drnkspwk]*m_drnkspwk if e(sample)
gen z3 = _b[_cons] + _b[3.physact]*1 + _b[age10 ]*m_age10 + _b[smoking]*m_smoking + _b[drnkspwk]*m_drnkspwk if e(sample) 
gen z4 = _b[_cons] + _b[4.physact]*1 + _b[age10 ]*m_age10 + _b[smoking]*m_smoking + _b[drnkspwk]*m_drnkspwk if e(sample) 
gen z5 = _b[_cons] + _b[5.physact]*1 + _b[age10 ]*m_age10 + _b[smoking]*m_smoking + _b[drnkspwk]*m_drnkspwk if e(sample) 
summ z1 z2 z3 z4 z5
* check that the differences between the adjusted marginal means 
* (relative to the reference activity level) equal the estimated regression coefficients
margins,  dydx(physact)

* test for heterogeneity
testparm i.physact
* test for linear trend
contrast q(1).physact, noeffects
* test for departure from linear trend
contrast q(2/4).physact, noeffects

* models using log of creatinine level as outcome and predictor
regress lncreat bmi age, eform("exp(beta)")
nlcom 100*(exp(_b[bmi])-1)
regress sbp lncreat age diabetes
nlcom _b[lncreat]*log(1.25)

log close