*** Factor analysis
factor  GSE1-GSE6,  ml
screeplot
rotate, oblimin

**** Calculating Alpha 
alpha GSE1-GSE6, item std  detail

*** Calculating bootstrapped confidence intervals for alpha 
bootstrap r(alpha), reps(1000): alpha GSE1-GSE6, item std detail
estat bootstrap, all

*** Calculating Omega after factor analysis (point estimate only, no confidence intervals)
*** Omega = (sum of factor loadings)^2 / (sum of uniquenesses + (sum of factor loadings)^2)
set more off
scalar drop _all
factor  GSE1-GSE6, factor(1) ml
rotate, oblimin 
matrix c = e(L)
matrix e = e(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 



* Calculating bootstrapped confidence intervals for Omega after EFA 

capture program drop omegaefa
program define omegaefa, rclass
version  12.0
set more off
preserve
syntax [varlist] [if] [in]
factor  `varlist', factor(1) ml
rotate, oblimin
matrix c = e(L)
matrix e = e(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 omega=c[1,1]/(c[1,1]+e[1,1])
restore
end

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