set memory 400M
set matsize 400

** use data from your own directory **
use "/Users/haley/Dropbox/Desktop/mlm09/data/zipcode_mortality.dta"

tab agecat
tab sex
tab y race, column
tab zip race

*** analysis assuming independence between individuals ***
xi:logit y race i.agecat sex
xi:logistic y race i.agecat sex
*** add interaction between age and sex ******
gen intagesex = agecat * sex
xi:logit y race i.agecat sex i.intagesex
xi:logistic y race i.agecat sex i.intagesex


*** random effects models ****

** install needed commands **
ssc describe gllamm
ssc install gllamm, replace /* to replace previous version */

** cannot use "xi:" with gllamm **
** type all the indicator variables for age category **
** keep the same age and interaction indicator variables ** 
** for consistent comparison with independence model results **

****************************
** random intercept model **
****************************
gllamm y race _Iagecat_2 _Iagecat_3 _Iagecat_4 _Iagecat_5 sex, i(zip) l(logit) f(binom) adapt
*** calculate average predicted probility ***
predict fit1,xb
gen phat1 = exp(fit1)/( 1 + exp(fit1) )
table race, c(mean phat1)

*** calculate predicted probability for each zip code ***
gllapred ebfit1, linpred
gen ebphat1 = exp(ebfit1)/( 1 + exp(ebfit1) )
table zip race, c(mean ebphat1)

gllamm y race _Iagecat_2 _Iagecat_3 _Iagecat_4 _Iagecat_5 sex _Iintagesex_1 _Iintagesex_2 _Iintagesex_3 _Iintagesex_5, i(zip) l(logit) f(binom) adapt

*** calculate average predicted probility ***
predict fit2, xb
gen phat2 = exp(fit2)/( 1 + exp(fit2) )
table race, c(mean phat2)

*** calculate predicted probability for each zip code ***
gllapred ebfit2, linpred
gen ebphat2 = exp(ebfit2)/( 1 + exp(ebfit2) )
table zip race, c(mean ebphat2)
 

**************************************
** random intercept and slope model **
**************************************

** add random slope to race coefficient **

** define equations used to specify random intercept and slopes **
gen byte one = 1  
eq cons : one  
eq race : race  

gllamm y race _Iagecat_2 _Iagecat_3 _Iagecat_4 _Iagecat_5 sex, i(zip) l(logit) f(binom) adapt nrf(2) eqs(race cons)

*** calculate average predicted probility ***
predict fit3,xb
gen phat3 = exp(fit3)/( 1 + exp(fit3) )
table race, c(mean phat3)

*** calculate predicted probability for each zip code ***
gllapred ebfit3, linpred
gen ebphat3 = exp(ebfit3)/( 1 + exp(ebfit3) )
table zip race, c(mean ebphat3)

** display zip code-specific predictions from 3 models together for easy comparison **
table zip race, c(mean ebphat1 mean ebphat2 mean ebphat3)