capture log close
log using myoc, replace 

set more off

** Data preparation
 * Original data
insheet using mi.dat, clear

 * List data
list q h l c count, noobs clean 

 * Prepare data
rename q y1                     /* y will be response            */
rename h y2
rename l y3
rename c y4
gen patt=_n                     /* identifier for level-2 units  */
gen wt2 = count                 /* counts are level-2 weights    */
reshape long y, i(patt) j(var)  /* stack all responses into y    */
tab var, gen(d)                 /* create dummy variables        */

 * list data for first response pattern
list patt var y d1 d2 d3 d4 wt2 if patt==1, clean


** Estimation in gllamm exploratory latent class model (as random coefficient model)
 * Equations 
eq d1: d1
eq d2: d2
eq d3: d3
eq d4: d4

* gllamm command
gllamm y, i(patt) ip(fn) nip(2) nrf(4) eqs(d1 d2 d3 d4) /*
  */ weight(wt) l(logit) f(binom) nocons


* Estimate equivalent two-class factor model
eq load: d1 d2 d3 d4
gllamm y d1 d2 d3 d4, i(patt) eqs(load) ip(f) nip(2) weight(wt) /*
  */ l(logit) f(binom) nocons iterate(30)


** Conditional response probabilities
 * Predictions based on two-class one-factor model (from exercise)
matrix m=e(zlc2)         /* locations e1 and e2 */
gen u01 = m[1,1]         /* class 1: e1 */
gen u11 = m[1,2]         /* class 2: e2 */
gllapred mu0, mu us(u0)  /* 1-specificity */
gllapred mu1, mu us(u1)  /* sensitivity   */

sort patt var
list var mu0 mu1 if patt==1, noobs clean


** Posterior probabilities using gllapred
gllapred prob, p


** Expected counts
replace wt2=1 if wt == 0  /* need l_j for patterns with wt2=0 */    
gllapred l, ll            /* log-likelihood contributions l_j */
gen expcount = 94*exp(l)


** Posterior probabilities and expected counts
format expcount %4.2f     /* format for list (2 dec. places) */
format prob2 %4.3f
sort patt
list patt count expcount prob2 if d1==1, clean

log close
exit