set more off
cd "~/Documents/teaching/c2021/biostat208/labs/lab4"
*do "~/Documents/teaching/c2021/biostat208/labs/lab4/lab4.do"

capture log close
log using lab4, replace


use lab4, clear

label list

* indicator for being relatively inactive
recode physact 2=1 3/5=0, gen(lessactive)
label variable lessactive "less active than peers"
label values lessactive noyes
tab physact lessactive

* three-level alcohol use variable
recode drnkspwk 0.2/4.9999 = 1 5/max = 2, gen(drinkamt)
label variable drinkamt "alcohol consumption"
label define drinkamt 0 "none" 1 "<5 drinks/week" 2 ">= 5 drinks/week"
label values drinkamt drinkamt
tab drnkspwk drinkamt

* natural log of triglyceride level
gen lntg = log(tgl)

* select for complete data so nested models run on same sample
foreach x in lncr bmi bmi age raceth educ drinkamt lessactive lntg {
	drop if missing(`x')
	}

* unadjusted model for crude association of BMI with log-creatinine levels
reg lncreat bmi, eform("exp(beta)")
* Relative increment in creatinine associated with 5-kg/m^2 increment in BMI
lincom bmi*5, eform
* Percent increment in creatinine associated with 5-kg/m^2 increment in BMI
nlcom 100*(exp(_b[bmi]*5)-1)

* model adjusting for confounders
regress lncreat bmi age i.raceth educyrs i.drinkamt lessactive, eform("exp(beta)")
* Relative increment in creatinine associated with 5-kg/m^2 increment in BMI
lincom bmi*5, eform
* Percent increment in creatinine associated with  5-kg/m^2 increment in BMI
nlcom 100*(exp(_b[bmi]*5)-1)

* Mediation
* Model 1:
reg lntg bmi age i.raceth educ i.drinkamt lessact
scalar link1 = _b[bmi]
* Percent increase in TG for a 5-kg/m^2 increase in BMI
nlcom 100*(exp(_b[bmi]*5)-1)

* Model 2:
reg lncr bmi age i.raceth educ i.drinkamt lessactive
* overall effect of BMI on creatinine levels
*  needed for evaluating mediation
scalar b_overall = _b[bmi]
* Overall BMI effect: percent increase in creatinine for a 5-kg/m^2 increase in BMI
nlcom 100*(exp(_b[bmi]*5)-1)

* Model 3:
reg lncr bmi age i.raceth educ i.drinkamt lessact lntg
scalar link2 = _b[lntg]
scalar b_direct = _b[bmi]
* Direct BMI effect: percent increase in creatinine for a 5-kg/m^2 increase in BMI
nlcom 100*(exp(_b[bmi]*5)-1)

* Relative increase in creatinine for a 25% increast in TG
local l125 = log(1.25)
lincom lntg*`l125', eform
* Percent increase in creatinine for a 25% increase in TG
nlcom 100*(exp(_b[lntg]*log(1.25))-1)

* indirect effect of BMI on creatinine via TG
*    on log-creatinine scale, for comparison with medeff results below
display link1*link2
*    percent increase in creatinine for a 5 kg/m^2 increase in BMI
display 100*(exp(5*link1*link2)-1)

* indirect effect via subtraction
display b_overall-b_direct
*    percent increase in creatinine for a 5-kg/m^2 increase in BMI
display 100*(exp(5*(b_overall-b_direct))-1)


* Percentage of adjusted BMI effect on log-creatinine explained by TG
scalar pe = (b_overall - b_direct)/b_overall*100
display round(pe, .1)

capture program drop mediate
program define mediate, rclass
	syntax varlist, outcome(varlist) mediator(varlist fv) covars(varlist fv)
	reg `outcome' `varlist' `covars'
	scalar b_overall = _b[`varlist']
	reg `outcome' `varlist' `covars' `mediator'
	scalar b_direct = _b[`varlist']
	return scalar pe = (b_overall - b_direct) / b_overall * 100
	end
* set seed for random number generator to make bootstrap results replicable
set seed 208
bootstrap pe=r(pe), reps(1000) notable nodots: mediate bmi, ///
    outcome(lncr) mediator(lntg) covars(age i.raceth educ i.drinkamt lessactive)
estat bootstrap, all

* Optional: estimates using medeff package
foreach x in raceth drinkamt {
	qui tab `x', gen(`x'_)
	drop `x'_1
	}
medeff (regress lntg bmi age raceth_* educ drinkamt_* lessact) /// // model 1
	(regress lncr bmi age raceth_* educ drinkamt_* lessact lntg), /// // model 3
	mediate(lntg) treat(bmi) sims(2000)

log close