---
title: "Bayesian Multilevel Factor Model (Part I)"
description: |
Comparing different approaches for fitting multilevel factor models with random intercepts
author: Mark Lai
date: "2024-07-31"
categories:
- Multilevel structural equation modeling
- Bayesian
---
\newcommand{\bv}[1]{\boldsymbol{\mathbf{#1}}}
\DeclareMathOperator{\Tr}{Tr}
In this note, I explore (mainly Bayesian) methods for fitting a multilevel factor model.
```{r}
library(lavaan)
library(cmdstanr)
register_knitr_engine()
library(blavaan)
future::plan("multicore", workers = 3)
```
## Model Equation
For response vector $\bv y_{ij}$ of length $p$ for person $i$ in cluster $j$, the multilevel factor model is
$$
\bv y_{ij} = \bv \nu + \bv \Lambda \bv \eta_{ij} + \bv \varepsilon^{b}_j + \bv \varepsilon^{w}_{ij}
$$
with
$$
\begin{aligned}
E(\bv \varepsilon^{b}_j) & = \bv 0 \\
E(\bv \varepsilon^{w}_{ij}) & = \bv 0 \\
E(\bv \eta_{ij}) & = \bv \alpha \\
V(\bv \varepsilon^{b}_j) & = \bv \Theta^b \\
V(\bv \varepsilon^{w}_{ij}) & = \bv \Theta^w \\
V(\bv \eta_{ij}) & = \bv \Psi^w + \bv \Psi^b
\end{aligned}
$$
::: {.callout-note}
## Identification Constraints
For each component in $\bv \eta$, we need one constraint on the mean structure and one constraint on the covariance structure. Conventionally, for the mean structure, we either set $E(\eta) = 0$ or $\nu = 0$ for one item; for the covariance structure, we either set $V(\eta) = 1$ or $\lambda = 1$ for one item.
:::
### Assumptions
Here we assume cross-level loading invariance (see [this paper](https://journals.sagepub.com/doi/abs/10.3102/1076998616646200) and [this paper](https://journals.sagepub.com/doi/pdf/10.1177/25152459231182319), for example), such that there is only one loading matrix $\bv \Lambda$ that applies to both levels, so that the latent variable has both a within-level and a between-level component:
$$
\bv \eta_{ij} = \bv \eta^b_j + \bv \eta^{w}_{ij}
$$
However, the loadings can still be randomly varying across clusters, which will be a discussion for Part II.
We also assume that $\bv \varepsilon^{b}_j$ is [i.i.d.](https://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables) across clusters and that $\bv \varepsilon^{w}_{ij}$ is i.i.d. across individuals.
::: {.callout-important}
## Conditional Independence
In Bayesian estimation we try to find conditional independence in the data so that the estimation can be simplified. Here, assuming multivariate normality of $\bv \varepsilon^{b}$ and $\bv \varepsilon^{w}$, we have independent observations when conditioned on the level-2 parameters:
$$
\bv y_{ij} \mid \bv \eta^b_j, \bv \varepsilon^{b}_j \sim N(\nu + \bv \varepsilon^{b}_j + \bv \Lambda \bv \alpha, \bv \Sigma^w),
$$
where $\bv \Sigma^w = \bv \Lambda \bv \Psi^w \bv \Lambda^\top + \bv \Theta^w$. Note that the within-cluster covariance is constant.
:::
## Simulate Data
Set up: One factor, Four items; two items with random intercepts
```{r}
# Design parameters
set.seed(123)
num_clus <- 100
clus_size <- 10
icc_eta <- .10
lambda <- c(.7, .9, .7, .8)
nu <- c(0, 0.3, -0.5, 0.2)
theta_lambda <- c(0, 0, 0, 0.08) # loading variances
theta_nu <- c(0, 0, 0.1, 0.05) # intercept variances
theta <- c(.51, .51, .40, .40) # uniqueness
# assume one-factor model holds for now
# Simulate latent variable
clus_id <- rep(1:num_clus, each = clus_size)
etab <- rnorm(num_clus, sd = sqrt(icc_eta / (1 - icc_eta)))
etaw <- rnorm(num_clus * clus_size)
eta <- etab[clus_id] + etaw
# confirm the simulated variable behaves as expected
# lmer(eta ~ (1 | clus_id))
# Simulate items for each cluster
num_items <- 4
y <- lapply(1:num_clus, FUN = \(j) {
yj <- nu + rnorm(nu, sd = sqrt(theta_nu)) +
tcrossprod(
lambda + rnorm(lambda, sd = sqrt(theta_lambda)),
etab[j] + etaw[clus_id == j]
) +
rnorm(num_items * clus_size, sd = sqrt(theta))
t(yj)
})
dat <- do.call(rbind, y) |>
as.data.frame() |>
setNames(paste0("y", 1:4)) |>
cbind(id = clus_id)
```
### Frequentist with `lavaan`
```{r}
mcfa_mod <- "
level: 1
f =~ l1 * y1 + l2 * y2 + l3 * y3 + l4 * y4
f ~~ 1 * f
level: 2
f =~ l1 * y1 + l2 * y2 + l3 * y3 + l4 * y4
"
mcfa_fit <- cfa(mcfa_mod, data = dat, cluster = "id", auto.fix.first = FALSE)
summary(mcfa_fit)
```
### Bayesian with `blavaan`
To my knowledge, the current version of `blavaan` does not support cross-level constraints.
```{r}
#| error: true
mcfa_bfit <- bcfa(mcfa_mod, data = dat, cluster = "id", auto.fix.first = FALSE)
```
It can fit a model without cross-level invariance:
```{r}
mcfa_uncon_mod <- "
level: 1
fw =~ y1 + y2 + y3 + y4
level: 2
fb =~ y1 + y2 + y3 + y4
"
mcfa_uncon_bfit <- bcfa(mcfa_uncon_mod, data = dat, cluster = "id",
bcontrol = list(cores = 3))
summary(mcfa_uncon_bfit)
```
It takes about `r mcfa_uncon_bfit@timing$Estimate` seconds to estimate the model.
## Using STAN
### Columnwise approach
Here's the code I used to fit multilevel factor model many years ago. It assumes that $\bv \Theta^w$ is diagonal, aka [*local independence*](https://en.wikipedia.org/wiki/Local_independence), so that I can treat each item as independent normal variable (after conditioning on the parameters). It does require sampling the latent variables: $\eta^b$, $\eta^w$, and $\varepsilon^b$.
::: {#lst-mcfa-col}
```{stan}
#| output.var: mcfa_col
#| file: "mcfa_col.stan"
```
STAN code for the columnwise approach
:::
```{r}
#| message: false
mcfa_col_fit <- mcfa_col$sample(
data = list(
N = nrow(dat),
p = 4,
Y = t(dat[1:4]),
J = length(unique(dat$id)),
jj = dat$id
),
chains = 3,
parallel_chains = 3,
iter_warmup = 500,
iter_sampling = 1000,
refresh = 500
)
```
```{r}
mcfa_col_fit$summary(c("lambda", "thetaw", "thetab", "psib")) |>
knitr::kable(digits = 2)
```
This is still relatively quick, but won't work when $\Theta^w$ is not diagonal.
### Marginal likelihood using sufficient statistics
From citations such as [this](https://link.springer.com/article/10.1007/bf02295842) and [this](https://www.tandfonline.com/doi/abs/10.1080/10705510903203466) and [this](https://www.mdpi.com/2624-8611/3/2/17), the likelihood function for cluster $j$ is
$$
\begin{aligned}
\ell_j & = n_j \log (2 \pi) + (n_j - 1) \log | \bv \Sigma^w_j | + \sum_{i = 1}^{n_j}(\bv y_{ij} - \bar{\bv y}_j)^\top {\bv \Sigma^w_j}^{-1} (\bv y_{ij} - \bar{\bv y}_j) \\
& \quad + \log | n_j \bv \Sigma^b + \bv \Sigma^w_j | + (\bar{\bv y}_j - \bar{\bv y})^\top (n_j \bv \Sigma^b + \bv \Sigma^w_j)^{-1} (\bar{\bv y} - \mu) + \\
& \quad (\bar{\bv y} - \mu)^\top {\bv \Sigma^w_j}^{-1} (\bv y_j - \bar{\bv y}_j)
\end{aligned}
$$ {#eq-llj}
This shows that the likelihood function has two multivariate normal components: that $(\bv y_j - \bar{\bv y}_j)$ is multivariate normal with mean 0 and covariance $\bv \Sigma^w_j$ with a sample size of $n_j - 1$, and that $(\bar{\bv y}_j - \bar{\bv y})$ is multivariate normal with mean 0 and covariance $n_j \bv \Sigma^b + \bv \Sigma^w_j$.
As $\bv a^\top \bv B \bv a$ is a scalar and can be written as $\Tr(\bv B \bv a \bv a^\top)$, the sufficient statistics in the above likelihood function is the cluster-specific sample cross-product matrix, $\bv S^w_j$ = $\sum_{i = 1}^{n_j} (\bv y_{ij} - \bar{\bv y}_j) (\bv y_{ij} - \bar{\bv y}_j)^\top$, and the sample cluster means, $\bar{\bv y}_j$.
::: {.callout-tip}
When the within-cluster covariance is assumed homogeneous, the likelihood function will depend only on the total within-sample cross-product matrix,
$$
\bv S^w_\text{pooled} = \sum_{j = 1}^J \bv S^w_j.
$$
:::
::: {.callout-tip}
When the design is balanced such that the cluster sizes are equal, on top of homogeneous of covariance condition, the likelihood function will depend only on the between-sample cross-product matrix,
$$
\bv S^b = \sum_{j = 1}^J (\bv y_j - \bar{\bv y}_j) (\bv y_j - \bar{\bv y}_j)^\top.
$$
:::
The marginal likelihood approach is also [used in `blavaan`](https://ecmerkle.github.io/blavaan/articles/multilevel.html). Below, I explore coding the model in STAN. Note that I ignore the mean structure since it is saturated in the factor model (as in single-level analysis), such that the last term in @eq-llj is 0.
#### Balanced design, pooled matrices
```{stan}
#| output.var: mcfa_suff_eq
#| file: "mcfa_suff_balanced.stan"
```
```{r}
#| message: false
# Prepare data for STAN
yc <- sweep(dat[1:4], MARGIN = 2, STATS = colMeans(dat[1:4])) # centered data
n <- table(dat$id)[1]
# Cross-product matrices
ssw <- crossprod(apply(yc, MARGIN = 2, FUN = \(x) x - ave(x, dat$id)))
ssb <- crossprod(apply(yc, MARGIN = 2, FUN = ave, dat$id))
# Run STAN
mcfa_suff_eq_fit <- mcfa_suff_eq$sample(
data = list(
p = 4,
J = num_clus,
n = n,
sw = ssw,
sb = ssb
),
chains = 3,
parallel_chains = 3,
iter_warmup = 500,
iter_sampling = 1000,
refresh = 500
)
```
```{r}
mcfa_suff_eq_fit$summary(c("lambda", "thetaw", "thetab", "psib")) |>
knitr::kable(digits = 2)
```
#### Using cluster-specific matrices
```{stan}
#| output.var: mcfa_suffj
#| file: "mcfa_suff2.stan"
```
```{r}
#| message: false
# Prepare data for STAN
nj <- table(dat$id)
# Cluster-specific cross-product matrices
sswj <- tapply(dat[1:4], INDEX = dat$id,
FUN = \(x) tcrossprod(t(x) - colMeans(x)))
ybarj <- tapply(dat[1:4], INDEX = dat$id, FUN = \(x) colMeans(x))
# ssbj <- tapply(yc, INDEX = dat$id,
# FUN = \(x) nrow(x) * tcrossprod(colMeans(x)))
# Run STAN
mcfa_suffj_fit <- mcfa_suffj$sample(
data = list(
p = 4,
J = num_clus,
n = nj,
sw = sswj,
ybar = ybarj
),
chains = 3,
parallel_chains = 3,
iter_warmup = 500,
iter_sampling = 1000,
refresh = 500
)
```
```{r}
mcfa_suffj_fit$summary(c("lambda", "thetaw", "thetab", "psib", "nu")) |>
knitr::kable(digits = 2)
```
Note that not even is the sampling faster, but the effective sample sizes are also higher, so the second approach is much more efficient.
```{r}
#| label: tbl-time
#| tbl-cap: "Sampling efficiency of different approaches"
data.frame(
model = c("blavaan", "columnwise", "sufficient statistics (balanced)", "sufficient statistics (cluster-specific)"),
ess = c(
blavInspect(mcfa_uncon_bfit, "draws") |>
posterior::summarise_draws("ess_bulk") |>
base::`[[`(2) |> base::`[`(16),
mcfa_col_fit$summary(variables = "psib", posterior::ess_bulk)[[2]],
mcfa_suff_eq_fit$summary(variables = "psib", posterior::ess_bulk)[[2]],
mcfa_suffj_fit$summary(variables = "psib", posterior::ess_bulk)[[2]]
),
time = c(
mcfa_uncon_bfit@timing$Estimate,
mcfa_col_fit$time()$total,
mcfa_suff_eq_fit$time()$total,
mcfa_suffj_fit$time()$total
)
) |>
knitr::kable(digits = 2, col.names = c("Model", "ESS (Bulk)", "Time (seconds)"))
```