capture log close
log using ichs, text replace
set more off
version 8.1

/*
Stata do-file for some of the analysis of Section 9.2 of
Skrondal, A. and Rabe-Hesketh, S (2004). Generalized
Latent Variable Modeling. Multilevel, Longitudinal
and Structural Equation Models. Chapman & Hall/CRC
*/

* read data
clear
infile id resp cons age xero cosine sine female height stunted time age1 season time2 using xerop.dat
list id time resp age xero cosine sine female height stunted in 1/20, clean

* Create indicator for last observation for each person
sort id time
qui by id: gen last=_n==_N

* GEE
xtgee resp age xero female cosine sine height stunted, i(id) corr(exch) l(logit) f(binom) robust

* save estimates for later
_estimates hold gee

* GEE with exponentiated coefficients
xtgee resp age xero female cosine sine height stunted, i(id) corr(exch) l(logit) f(binom) robust eform

* Random intercept model
gllamm resp age xero female cosine sine height stunted, i(id) l(logit) f(binom) nip(12) adapt

* Exponentiated coefficients
gllamm, eform

* Empirical Bayes predictions (to be stored in ebm1)
gllapred eb, u

* Simulate three datasets and obtain empirical Bayes predictions for Figure 9.2 
set seed 1322411
matrix aa=e(b)
gllasim y, from(aa)
gllamm y age xero female cosine sine height stunted, i(id) l(logit) f(binom) nip(12) from(aa) copy adapt
gllapred SIM1, u
drop y
gllasim y, from(aa)
gllamm y age xero female cosine sine height stunted, i(id) l(logit) f(binom)  nip(12) from(aa) copy adapt
gllapred SIM2, u
drop y
gllasim y, from(aa)
gllamm y age xero female cosine sine height stunted, i(id) l(logit) f(binom) nip(12) from(aa) copy adapt
gllapred SIM3, u
drop y

* Make Figure 9.2 and store as ichs2.gph

graph box ebm1 SIM1m1 SIM2m1 SIM3m1 if last==1, t2(" ") saving(ichs2, replace) /*
  */ legend(lab(1 "Data") lab(2 "Sim 1") lab(3 "Sim 2") lab(4 "Sim 3") ) /*
  */ marker(1, m(oh)) marker(2, m(oh)) marker(3, m(oh)) marker(4, m(oh)) lintensity(255)

* translate ichs2.gph ichs2.eps, translator(Graph2eps) replace /* to get postscript file */

* Make Figure 9.3 and store as ichs3.gph

sort ebm1
gen rank = sum(last)
serrbar ebm1 ebs1 rank if last==1&mod(rank,6)==5, xlab(5 100 200 275)  /*
  */ ytitle("Predicted random intercept") saving(ichs3, replace)

* translate ichs3.gph ichs3.eps, translator(Graph2eps) replace /* to get postscript file */

* Manipulate data, etc. for Figure 9.1

* Create small dataset with covariate values for which we want predictions
keep in 1/83
replace age = -32+_n-1
replace cosine=-1
replace sine= 0
replace xero=0
replace female=1
summ height
replace height=r(mean)
replace stunted=0
gen y=1

* Obtain marginal probabilities
gllapred mu_m, mu marg from(aa)
gen u01=0
gen u8p1 = 0.8
gen u8m1 = -0.8
gllapred mu_0, mu us(u0)
gllapred mu_8p, mu us(u8p)
gllapred mu_8m, mu us(u8m)
_estimates unhold gee
predict mu_g, mu
label variable mu_g "GEE"
label variable mu_0 "u=0"
label variable mu_8p "u=+0.8"
label variable mu_8m "u=-0.8"
label variable mu_m "pop. average"

sort age

* Make Figure 9.1 (ichs.gph)

twoway (line mu_m age, clpatt(solid)) /*
  */ (scatter mu_g age if mod(age,2)==1, ms(oh) mc(black)) (line mu_0 age, clpatt(dash)) /*
  */ (line mu_8p age, clpatt(dash)) (line mu_8m age, clpatt(dash)), /*
  */ legend(order(1 2 3) label(3 "conditional") position(1) ring(0) ) /*\
  */ ytitle("Probability of infection") /*
  */ xtitle("Age in months") saving(ichs, replace) 


* translate ichs.gph ichs.eps, translator(Graph2eps) replace /* to get postscript file */

log close
exit