/*** DAY 5 CODE ***/

*** 6.1 Single level models for dichotomous responses

** 6.1.1 Generalized linear model formulation
use http://www.stata-press.com/data/mlmus2/womenlf, clear
recode workstat 2=1
logit workstat husbinc chilpres 
logit workstat husbinc chilpres, or

predict prob, pr

twoway (line prob husbinc if chilpres==0, sort)         ///
 (line prob husbinc if chilpres==1, sort lpatt(dash)), ///
 legend(order(1 "No child" 2 "Child"))                  ///
 xtitle("Husband's income / $1000")                     ///
 ytitle("Probability that wife works")

twoway (function y=invlogit(_b[husbinc]*x+_b[_cons]), range(-100 100))    ///
 (function y=invlogit(_b[husbinc]*x+_b[chilpres]+_b[_cons]),              ///
 range(-100 100) lpatt(dash)),                                           ///
 xtitle("Husband's income / $1000") ytitle("Probability that wife works") ///
 legend(order(1 "No child" 2 "Child")) xline(1) xline(45)

glm workstat husbinc chilpres, link(logit) family(binom)

** 6.1.2 Latent-respons formulation
probit workstat husbinc chilpres

twoway (function y=invlogit(1.3358-0.0423*x)), range(-100 100))           ///
 (function y=norm(0.7982-0.0242*x), range(-100 100) lpatt(dash)),        ///
 xtitle("Husband's income / $1000") ytitle("Probability that wife works") ///
 legend(order(1 "Logit link" 2 "Probit link")) xline(1) xline(45)


*** 6.3 Longitudinal data structure

use http://www.stata-press.com/data/mlmus2/toenail, clear

xtdescribe if outcome<., i(patient) t(visit)


*** 6.4 Population-averaged or marginal probabilities

label define tr 0 "Itraconazole" 1 "Terbinafine"
label values treatment tr

graph bar (mean) proportion= outcome, over(visit) by(treatment) ///
   ytitle(Proportion with onycholysis)

egen prop = mean(outcome), by(treatment visit)
egen mn_month = mean(month), by(treatment visit)

twoway line prop mn_month, by(treatment) sort ///
  xtitle(Time in months) ytitle(Probability of onycholysis)

generate trt_month = treatment*month
logit outcome treatment month trt_month, or
predict prob, p

twoway (line prop mn_month, sort) (line prob month, sort lpatt(dash) ),          ///
   by(treatment) legend(order(1 "Observed proportions" 2 "Fitted probabilities")) ///
   xtitle(Time in months) ytitle(Probability of onycholysis)


*** 6.7 Estimation of logistic random intercept model

** 6.7.1 Using xtlogit
xtlogit outcome treatment month trt_month, i(patient) quad(30)
xtlogit, or

** 6.7.2 Using xtmelogit
xtmelogit outcome treatment month trt_month || patient: , intpoints(30)
estimates store xtmelogit

** 6.7.3 Using gllamm
gllamm outcome treatment month trt_month, i(patient) link(logit) family(binom) ///
 nip(30) adapt
estimates store gllamm
gllamm, eform


*** 6.8 Inference for random-intercept logistic regression 

lincom treatment + trt_month*20
lincom treatment + trt_month*20, or
lincom -(treatment + trt_month*20), or


*** 6.9 Subject-specific vs. population averaged relationship

display (invlogit(1) + invlogit(2))/2
display invlogit((1+2)/2)


*** 6.10 Measures of dependence and heterogeneity

** 6.10.2 Median odds ratio

display exp(sqrt(2*16.08)*invnorm(3/4))


*** 6.11 Maximum likelihood estimation

** 6.11.2 Some speed considerations

xtmelogit outcome treatment month trt_month || patient: , intpoints(1)


*** 6.12 Assigning values to random effects

** 6.12.1 Maximum likelihood estimation

estimates restore xtmelogit
predict offset, xb

statsby mlest=_b[_cons], by(patient) saving(ml): logit outcome, offset(offset)

sort patient
merge patient using ml
drop _merge