Higher-order IRT model in brms

Is it possible to estimate a higher-order item response theory (IRT) model in brms? Here are some example references: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.990.199&rep=rep1&type=pdf ; and https://www.frontiersin.org/articles/10.3389/fpsyg.2019.01328/full.

Below, I provide a reproducible example of a multidimensional IRT model, per the discussion here (Multidimensional IRT in brms), that I’d like to extend to estimate a higher-order latent factor. In the multidimensional model example below, there are three dimensions (A, B, and C). I’d like to extend the model to allow the three dimensions to load onto a higher-order latent factor (g), where g reflects the common variance among the dimensions. That is, g influences each of the dimensions, which in turn, influence people’s responses on the items (similar to what is depicted here: https://upload.wikimedia.org/wikipedia/commons/5/50/Carroll_three_stratum_model_of_human_Intelligence.png). Basically, I’d like to get estimates of people’s levels (theta) on each dimension in addition to their level on the overall factor. Is this possible in brms? If so, how can I modify the model below to achieve this?

Thanks in advance!

library("brms")
library("mirt")
library("mvtnorm")
library("tidyverse")

numberItems <- 30
numberItemLevels <- 2
numberDimensions <- 3
sampleSize <- 1000

#Simulate Data
set.seed(1)

itemsPerDimension <- ceiling(numberItems/numberDimensions)

variances <- rep(10, numberDimensions)

corrs <- matrix(.5, nrow = numberDimensions, ncol = numberDimensions)
diag(corrs) <- 1

covar <- outer(sqrt(variances), sqrt(variances))*corrs

d <- matrix(rnorm(numberItems * numberItemLevels), nrow = numberItems)
d <- t(apply(d, 1, sort, decreasing = TRUE))

a <- matrix(NA, nrow = numberItems, ncol = numberDimensions)
for(i in 1:numberDimensions){
  a[(i*itemsPerDimension - itemsPerDimension + 1):(i*itemsPerDimension), i] <- 1
}

simulatedData <- simdata(a, d, sigma = covar, N = sampleSize, itemtype = "graded")

#Wide to Long form
dat <- data.frame(simulatedData)
dat$id <- as.factor(as.numeric(row.names(dat)))

dat_long <- pivot_longer(dat, -id, names_to = c("item"), values_to = "response")

#Determine item dimensions
dat_long$itemNumber <- parse_number(dat_long$item)
dat_long$dimension <- factor(LETTERS[ceiling(dat_long$itemNumber/itemsPerDimension)])

#Change data types
dat_long$item <- factor(dat_long$item)
dat_long$response <- factor(dat_long$response, levels = c(0,1,2), ordered = TRUE)

#Model formula
modelFormula <- bf(
  response ~ 1 + dimension + (1 | item) + (1 + dimension | id),
  disc ~ 1 + dimension (1 | item)	
)

#Fit the model
modelFit <- brm(
  formula = modelFormula,
  data = dat_long,
  family = brmsfamily("cumulative", "logit"),
  seed = 1234
)

As far as I know one cannot estimate higher level “factors” with brms “out of the box”.
But @paul.buerkner knows this better.
If you are familiar with writing models in Stan, you can maybe extend a brms model by using brms` stanvars functionality to add data, parameters and model specifications to a brms model.

See here for a related question: Example hierarchical models with 3+ levels of hierarchy?