#Read  data Mplus example  bi-factor EFA
ex5.29 <- read.table("http://statmodel.com/usersguide/chap5/ex5.29.dat")
names(ex5.29) = paste0("y",1:10)
#install packages
library(lavaan)
library(lavaanPlot)
library(GPArotation)
library(semTools)

#Define EFA block
model.BESEM <- '
  efa("efa1")*FG +
  efa("efa1")*F1 +
  efa("efa1")*F2 =~ y1+y2+y3+y4+y5+y6+y7+y8+y9+y10
  '
#fit.BESEM
fit.BESEM <- cfa(
  model = model.BESEM,
  data = ex5.29,
  estimator = "ML",
  rotation = "bigeomin",
  rotation.args = list(
    rstarts = 30,
    algorithm = "gpa",
    orthogonal = TRUE,
    std.ov = TRUE,
    geomin.epsilon = 0.0001
  ),
)
  summary(fit.BESEM, standardized = TRUE)
  #BifactorESEM graph 
  lavaanPlot(name = "BifactorESEM",
             model = fit.BESEM)