sjlog using diet, replace
set more off

/* cme command */
use diet_w, clear

gen lfiber1 = ln(fiber1)-2.8
gen lfiber2 = ln(fiber2)-2.8

cme chd bus (lfib: lfiber1 lfiber2), l(logit) f(binom) nolog
cme, eform ind total


/* sensitivity analysis */
gen beta_u = .
gen beta_1 = .
gen variance = (_n-1)/200 in 1/11
glm chd bus lfiber1, l(logit) f(binom)
replace beta_u = [chd]lfiber1 in 1
replace beta_1 = [chd]bus in 1
forvalues i=2/11{
	local var=(`i'-1)/200
	qui cme chd bus (lfib: lfiber1), l(logit) f(binom) mevar(`var')
	replace beta_u = [chd]lfib in `i'
	replace beta_1 = [chd]bus in `i'
}

set scheme sj

twoway connect beta_1 variance, ytitle(Log odds ratio for bus) /*
  */ xtitle(Measurement error variance)
graph export sensit1.eps, replace


twoway connect beta_u variance, ytitle(Log odds ratio for true log-fiber) /*
  */ xtitle(Measurement error variance)
graph export sensit.eps, replace


cme chd (lfib: lfiber1 lfiber2), l(logit) f(binom) nolog tcovmod(bus)

/* using gllamm */
cme chd bus (lfib: lfiber1 lfiber2), l(logit) f(binom) nolog commands


*------------ commands created by cme ----------------------
* starting values
matrix startv = ( -.1474, -1.853, -1.458, -1.54, .3555, -.1237, .06659)
gen _id = _n

* collapse data to make gllamm faster
gen _one = 1
collapse (sum) _wt2 = _one, by(chd lfiber1 lfiber2 _id bus bus)

* give response variable and replicate measurements same prefix
rename chd _r1
rename lfiber1 _r2
rename lfiber2 _r3

* reshape data to long
reshape long _r, i(_id) j(_var)

* create dummy variables and interactions
gen byte cons = 1
gen byte _d1 = _var == 1
gen byte _dmeas = 1-_d1
gen _type = _d1 + 2*_dmeas /* response type */
gen _bus_d1 = bus*_d1

* define equations
eq load: _dmeas _d1
eq f1: bus cons
 
* call gllamm
gllamm _r  _bus_d1 _d1, /* 
 */ i(_id) nocons eqs(load) link(logit ident) family(binom gauss) /*
 */ lv(_type) fv(_type)  geqs(f1)  from(startv) copy adapt  /*
 */  weightf(_wt)        nolog 

*------------------- end of commands  ------------------------------


/* NPMLE (direct and indirect effects)  */

eq load: _dmeas _d1
eq f1: bus cons

gllamm _r _bus_d1 _d1, /* 
 */ i(_id) nocons eqs(load) link(logit ident) family(binom gauss) /*
 */ fv(_type) lv(_type) geqs(f1) nip(3) ip(f)

global HG_error=0
local n = 4
while $HG_error==0{
	matrix a=e(b)
	local k = e(k)
	local ll = e(ll)
	gllamm _r  _bus_d1 _d1, /* 
	 */ i(_id) nocons eqs(load) link(logit ident) family(binom gauss) /*
	 */ fv(_type) lv(_type) geqs(f1) nip(`n') ip(f) lf0(`k' `ll') from(a) gateaux(-2 2 1000)
	local n= `n'+1
}

* save locations and log probabilities as variables
matrix locs=e(zlc2)'
matrix lp = e(zps2)'
svmat locs 
svmat lp

* calculate probabilities
gen p=exp(lp1)

* plot masses
twoway (dropline p locs1), xtitle(Location) ytitle(Probability)
graph export mass.eps, replace

* set end-points for cumulative probabilities
replace locs1 = -1 in 9
replace p = 0 in 9
replace locs1 = 1.1 in 10
replace p = 0 in 10

* calculate cumulative probabilities
sort locs1
gen cump = sum(p)

* plot cumulative probabilties
twoway (line cump locs1, connect(stairstep))              /*
  */   (function y=norm(x/sqrt(.07019772)), range(-1 1)), /*
  */   xtitle(Location) ytitle(Cumulative probability)    /*
  */   legend(order(1 "NPMLE" 2 "Normal"))
graph export cum.eps, replace


/******* congeneric measurement model ********************/
expand 2 if _var==3
qui tab _var, gen(d)
bysort _id _var: gen d4 = _n==2
replace d3=0 if d4==1
eq load: d2 d3 d4 d1
eq het: d2 d3 d4
eq f1: bus

matrix a=J(1,12,0)
qui gllamm _r d1 _dmeas d3 d4, /* 
 */ i(_id) nocons eqs(load) link(logit ident) family(binom gauss) /*
 */ fv(_type) lv(_type) geqs(f1) s(het) from(a) copy eval
matrix a=e(b)
matrix list a
* Intercept for CHD
matrix a[1,1]= - 2
* mu
matrix a[1,2] = 3  /* was or 0? */
* Bias
matrix a[1,3]=1       /* measure 2 */
matrix a[1,4]=-1      /* measure 3 */
* log-sd
matrix a[1,5]=-2     /* measure 1 */
matrix a[1,6]=-2     /* measure 2 */
matrix a[1,7]=-1     /* measure 3 */
* factor loadings
matrix a[1,8]=1.5 /* measure 2 */
matrix a[1,9]=2   /* measure 3 */
matrix a[1,10]=-2
* true covariate sd
matrix a[1,11]=0.3
* effect of bus on true covariate
matrix a[1,12]= -0.2
matrix list a

set seed 131123
gllasim r, fsample from(a)
gllamm r _d1 _dmeas d3 d4, /* 
 */ i(_id) nocons eqs(load) link(logit ident) family(binom gauss) /*
 */ fv(_type) lv(_type) geqs(f1) s(het) adapt nip(12)

eq load: _dmeas _d1
gllamm r _d1 _dmeas d3 d4, i(_id) nocons eqs(load) link(logit ident) /*
  */ family(binom gauss) fv(_type) lv(_type) geqs(f1) adapt nip(12)

sjlog close, replace
exit