set more off
* gllamm flavor: estimate or use?
args flavor
* default is to use the existing results

* checking the results against -gllamm-
use hs-cfa

* estimation by -confa-
confa (vis: x1 x2 x3) (text: x4 x5 x6) (math: x7 x8 x9), from(iv)
est store confa_l
est save confa_l, replace
mat b_confa_l = e(b)
mat V_confa_l = e(V)

confa (vis: x1 x2 x3) (text: x4 x5 x6) (math: x7 x8 x9), from(b_confa_l) robust search(off)
est store confa_r
est save confa_r, replace
mat V_confa_r = e(V)

if "`flavor'" == "estimate" {

     * setting the data up for -gllamm-
     reshape long x, i(id) j(item 1-9) /* xi(sex age* school grade) */
     tab item, gen( i )
     eq f1 : i1 i2 i3
     eq f2 : i4 i5 i6
     eq f3 : i7 i8 i9
     eq theta : i1 i2 i3 i4 i5 i6 i7 i8 i9
     reshape long x, i(id) j(item 1-9)
     tab item, gen( i )

     * starting values: all ones
     mat J30 = J(1,30,1)
     * gllamm with the loadings standardization
     gllamm x i1-i9, eq(f1 f2 f3) nrf(3) nip(4 4 4) i(id) s(theta) ///
          nocons trace from( J30 ) copy dots adapt difficult
     est store gllamm_nip4_l
     est save gllamm_nip4_l, replace
     mat b_gllamm_nip4_l = e(b)

     * increase the number of integration points
     gllamm x i1-i9, eq(f1 f2 f3) nrf(3) nip(12 12 12) i(id) s(theta) ///
          nocons trace from( b_gllamm_nip4_l ) copy dots
     est store gllamm_nip12_l
     est save gllamm_nip12_l, replace

     * -robust- option
     gllamm x i1-i9, eq(f1 f2 f3) nrf(3) nip(12 12 12) i(id) s(theta) ///
          nocons trace from( b_gllamm_nip4_l ) copy dots robust
     est store gllamm_nip12_r
     est save gllamm_nip12_r, replace

}

est use gllamm_nip12_l
mat b_gllamm_nip12_l = e(b)
mat V_gllamm_nip12_l = e(V)
* store some results and their transformations: factor variances
mat chol_g = e(chol)
* and error variances
mat th_g = 0
mat th_V = 0
forvalues k=1/9 {
  nlcom (theta_`k': exp( [lns1]i`k' ) )
  mat th_g = th_g, r(b)
  mat th_V = th_V , r(V)
}
mat th_g = th_g[1,2..10]
mat th_V = th_V[1,2..10]

est use gllamm_nip12_r
mat V_gllamm_nip12_r = e(Vr)
mat th_Vr = 0
forvalues k=1/9 {
  nlcom (theta_`k': exp( [lns1]i`k' ) )
  mat th_Vr = th_Vr , r(V)
}
mat th_Vr = th_Vr[1,2..10]

mata
// POINT ESTIMATES
CONFA = st_matrix( "b_confa_l" )
GLLAMM = st_matrix( "b_gllamm_nip12_l" )

// means; scm means "Select from Confa, Means", etc.
scm = J(1,9,0), J(1,24,0)
sgm = J(1,9,0), J(1,21,0)
mreldif( select(CONFA,scm), select(GLLAMM,sgm))

// lambdas
scl = J(1,10, 0) , 1, 1, 0, 1, 1, 0, 1, 1, J(1, 15, 0)
sgl = J(1,18,0) , 1, 1, 0, 1, 1, 0, 1, 1, J(1, 4, 0)
mreldif( select(CONFA,scl), select(GLLAMM,sgl))

// factor variances
scp = J(1,18,0), J(1,6,1), J(1,9,0)
chol_g = st_matrix("chol_g")
Phi_g = chol_g * chol_g'
vechPhi_g = Phi_g[1,1], Phi_g[2,2], Phi_g[2,1], Phi_g[3,3], Phi_g[3,2], Phi_g[3,1]
mreldif( select(CONFA,scp), vechPhi_g )

// error variances
sct = J(1,24,0), J(1,9,1)
th_g = st_matrix("th_g")
mreldif( select(CONFA,sct), th_g )

// STANDARD ERRORS
V_CONFA = st_matrix("V_confa_l")
V_GLLAMM = st_matrix("V_gllamm_nip12_l")

V_CONFA_r = st_matrix("V_confa_r")
V_GLLAMM_r = st_matrix("V_gllamm_nip12_r")

// means
mreldif( sqrt(diagonal(V_CONFA[1..9,1..9])), sqrt(diagonal(V_GLLAMM[1..9, 1..9])) )
mreldif( sqrt(diagonal(V_CONFA_r[1..9,1..9])), sqrt(diagonal(V_GLLAMM_r[1..9, 1..9])) )

// loadings
VLG = select( select( V_GLLAMM, sgl), sgl' )
VLC = select( select( V_CONFA, scl), scl' )
mreldif( sqrt(diagonal( VLG )), sqrt(diagonal( VLC )) )

VLG_r = select( select( V_GLLAMM_r, sgl), sgl' )
VLC_r = select( select( V_CONFA_r, scl), scl' )
mreldif( sqrt(diagonal( VLG_r )), sqrt(diagonal( VLC_r )) )

// factor variances, hand-copied from gllamm
VPG = (.15017148, .11083428, .07495402, .0919132, .04863624, .05474253)
VPC = sqrt(diagonal(select( select( V_CONFA, scp), scp' ) ))
mreldif( VPG, VPC' )

VPG_r = (.1821103, .11926282, .09050493, .10692003, .05493861, .05855645)
VPC_r = sqrt(diagonal(select( select( V_CONFA_r, scp), scp' ) ))
mreldif( VPG_r, VPC_r' )

// error variances
VTC = sqrt(diagonal(select( select( V_CONFA, sct), sct' ) ))
VTG = st_matrix("th_V")
mreldif( VTG, VTC' )

VTC_r = sqrt(diagonal(select( select( V_CONFA_r, sct), sct' ) ))
VTG_r = st_matrix("th_Vr")
mreldif( VTG_r, VTC_r' )

end

exit