***Calculating Omega after CFA with a single factor model using SEM
set more off
scalar drop _all
sem (L -> GSE1-GSE6, ),  latent(L) variance(L@1) standardized nocapslatent
estat framework
matrix c = r(Gamma)
matrix e = vecdiag(r(Psi))'
* sum(loadings)^2 
matrix one = J(rowsof(c),1,1)
matrix c = one'*c*c'*one

* sum(error variances)
matrix e = J(1,rowsof(e),1)*e

matrix list c
matrix list e
scalar omega= c[1,1]/(c[1,1]+e[1,1])
scalar list 


** Omega bootstrapped confidence intervals after CFA with a single factor model 

capture program drop omegasem
program define omegasem, rclass
version  12.0
set more off
preserve
syntax [varlist] [if] [in]
sem (L -> `varlist'),  latent(L) variance(L@1) standardized nocapslatent
estat framework
matrix c = r(Gamma)
matrix e = vecdiag(r(Psi))'

* sum(loadings)^2 
matrix one = J(rowsof(c),1,1)
matrix c = one'*c*c'*one

* sum(error variances)
matrix e = J(1,rowsof(e),1)*e

matrix list c
matrix list e
return scalar omegah=c[1,1]/(c[1,1]+e[1,1])
restore
end

bootstrap omega=r(omegah), reps(1000): omegasem GSE1-GSE6
estat bootstrap, all