```{r}
# install.packages(c('lavaan', 'psych', 'semPlot'))
library(lavaan)
library(psych)
library(semPlot)
```

```{r}
n <- 11914

# 1. Specify correlation matrix
cor.mat <- matrix(
  c(
    1, .722, .801, .689, .524, .452, .637, .695, .695, .760, .800,
    .722, 1, .708, .672, .627, .515, .533, .827, .684, .658, .727,
    .801, .708, 1, .803, .617, .550, .529, .670, .593, .684, .981,
    .689, .672, .803, 1, .608, .561, .423, .637, .521, .573, .901,
    .524, .627, .617, .608, 1, .701, .306, .617, .408, .421, .643,
    .452, .515, .550, .561, .701, 1, .225, .520, .336, .342, .578,
    .637, .533, .529, .423, .306, .225, 1, .415, .741, .745, .519,
    .695, .827, .670, .637, .617, .520, .415, 1, .600, .585, .688,
    .695, .684, .593, .521, .408, .336, .741, .600, 1, .743, .596,
    .760, .658, .684, .573, .421, .342, .745, .585, .743, 1, .679,
    .800, .727, .981, .901, .643, .578, .519, .688, .596, .679, 1
    ),
  nrow=11,
  byrow=T
)
colnames(cor.mat) <- rownames(cor.mat) <- c('GS', 'AR', 'WK', 'PC', 'NO', 'CS', 'AS', 'MK', 'MC', 'EI', 'VE')

# Subtests that compose the ASVAB and AFQT. VE (verbal expression) is just a composite of WK and PC), so it should not be included.
asvab <- c('GS', 'AR', 'WK', 'PC', 'NO', 'CS', 'AS', 'MK', 'MC', 'EI')
afqt <- c('AR', 'WK', 'PC', 'MK')
```

```{r}
fa.parallel(x=cor.mat[asvab, asvab], n.obs=n, fm='ml', n.iter=200)
```

```{r}
hier_model <- '
  # first‐order factors
    math   =~ AR + MK
    verbal =~ WK + PC
    tech   =~ AS + MC + EI + GS
    speed  =~ NO + CS

  # second‐order g
    g =~ math + verbal + tech + speed
'

fit <- cfa(
  hier_model,
  sample.cov  = cor.mat[asvab, asvab],
  sample.nobs = n,
  estimator = "ML"
)

summary(fit, standardized=T, rsquare=T)
stdf <- inspect(fit, "std")$lambda
print(stdf)
lavInspect(fit, 'cor.lv')
```

```{r}
semPaths(fit, what="std", edge.color="black", edge.label.cex=.9, esize=1, nDigits=2, sizeLat=7)
```
```{r}
g_efa <- omega(
  cor.mat[asvab, asvab],
  n.iter = 1,
  nfactors = 4,
  digits   = 3,
  n.obs    = n,
  sl       = T
) # sl=TRUE returns Schmid‐Leiman

print(g_efa, digits=3)
```

```{r}
hier_model <- '
  g =~ AR + MK + WK + PC + AS + MC + EI + NO + CS + GS
  math =~ AR + MK
  verbal =~ WK + PC
  tech =~ AS + MC + EI + GS
  speed =~ NO + CS
'

fit <- cfa(
  hier_model,
  sample.cov  = cor.mat[asvab, asvab],
  sample.nobs = n,
  estimator = 'ML',
  orthogonal=T,
  std.lv=T
)

g_cfa <- omegaFromSem(fit)
print(g_cfa, digits=3)
```

```{r}
fa.parallel(x=cor.mat[afqt, afqt], n.obs=n, fm='ml', n.iter=200)
```
```{r}
g_efa <- omega(
  cor.mat[afqt, afqt],
  nfactors  = 2,
  digits    = 3,
  n.obs     = n,
  sl        = T,
  estimator = 'ML'
)

print(g_cfa, digits=3)
```

```{r}
hier_model <- '
  # first‐order factors
    math   =~ AR + MK
    verbal =~ WK + PC
  # second‐order g
    g =~ a*math + a*verbal
'

fit <- cfa(
  hier_model,
  sample.cov  = cor.mat[afqt, afqt],
  sample.nobs = n,
  estimator = 'ML'
)

summary(fit, standardized=T, rsquare=T)

stdf <- inspect(fit, 'std')$lambda
print(stdf)
```

```{r}
semPaths(fit, what='std', edge.color='black', edge.label.cex=1, esize=1, nDigits=3, sizeLat=7)
```