higherorder.R

library("lavaan")
library("psychonetrics")
library("dplyr")

# Generate data:
mod <- '
F1 =~ X1 + X2 + X3
F2 =~ X4 + X5 + X6
F3 =~ X7 + X8 + X9
G =~ F1 + F2 + F3
'
dat <- simulateData(mod)

# Factor loadings matrix:
Lambda <- cbind(
  simplestructure(rep(1:3,each=3)),
  0
)

# beta matrix:
Beta <- matrix(0,4,4)
Beta[1:3,4] <- 1


# Make model:
mod <- lvm(dat, 
            lambda = Lambda, 
            beta = Beta,
            sigma_zeta = "diag",
            identification = "variance")

# Run model:
mod <- mod %>% runmodel

# Parameters:
mod %>% parameters

# Fit: 
mod %>% fit