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

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


* Read data

insheet using kenkel.dat, clear
list in 1/10, clean

* turn counts into integer values

replace drinks=round(drinks,1)


/* easy way to estimate the model, requires gllamm wrapper ssm (ssc install ssm) */

ssm drinks advice black hieduc, s(advice = black hieduc hlthins regmed heart)
  adapt q(16) family(poiss) link(log)

/* harder way to estimate model using gllamm directly */

* Collapse data and generate frequency weight variable wt2 to speed up estimation

disp _N
gen one=1
collapse (sum) wt2=one, by(black hlthins regmed heart hieduc drinks advice)
disp _N

gen id=_n
gen cons=1
list in 1/10, clean

* Poisson model for drinking

gllamm drinks advice cons hieduc black, i(id) weight(wt) family(poisson) link(log) nocons init 


* Overdispersed Poisson

gllamm drinks advice cons hieduc black, i(id) weight(wt) family(poisson) link(log) /*
    */ nocons adapt nip(10) 


* Probit model for advice

gllamm advice cons hieduc black hlthins regmed heart, i(id) weight(wt) family(binomial) /*
    */ link(probit) nocons init


* Prepare data for endogeneous treatment model 

* stack drinking and advice into single variable resp and create
* type = 1 for drinking and type = 2 for advice

rename drinks resp1 
gen resp2 = advice
reshape long resp, i(id) j(type)

* create dummies d1 for type=1 and d2 for type=2
tab type, gen(d)

sort id type
list in 1/10, clean

* create interactions between d1 and covariates in drining model

gen d1_advice = d1*advice
gen d1_hieduc = d1*hieduc
gen d1_black = d1*black

* create interactions between d2 and covariates in advice model
* (use foreach to save typing)

foreach var in hieduc black hlthins regmed heart {
    gen d2_`var' = d2*`var'
}

* Endogenous treatment: 

eq fac: d2 d1
gllamm resp d1_advice d1 d1_hieduc d1_black d2 d2_hieduc d2_black d2_hlthins /*
    */ d2_regmed d2_heart, nocons i(id) /*
    */ weight(wt) family(poisson binom) link(log probit) fv(type) lv(type) /*
    */ eq(fac) adapt nip(15) trace 


log close
exit