Code
library(MplusAutomation)
library(posterior)
library(cmdstanr)
register_knitr_engine()Bayesian Multilevel Factor Model (Part III)
Multilevel structural equation modeling
Bayesian
Fitting multilevel factor models with random factor loadings
Note
The code in this blog post was originally drafted in 2024, but then I got distracted with some other projects and personal matters. I’m finally able to get this done.
In Part I, I introduced the multilevel factor model and how to fit it using blavaan and Stan. In Part II, I explored the issue of heterogeneous latent variances across clusters and how to model them. In this Part III, I will further generalize the model to allow for random factor loadings, meaning that the relationship between the latent factor and the observed items can vary across clusters.
Simulate Data
Same set up: One-factor four items; heterogeneous latent variance; two items with random intercepts, one item with random loading
Code
# Design parameters
set.seed(1234)
num_clus <- 500
clus_size <- 10
icc_eta <- .10
lambda <- c(.7, .9, .7, .8)
nu <- c(0, 0.3, -0.5, 0.2)
range_dpsiw <- .80
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)))
psiw <- 1 + seq(-range_dpsiw, range_dpsiw, length.out = num_clus)
# Arrange etaw and psiw so that they're negatively correlated
etab <- sort(etab)[c(
sample((num_clus / 2 + 1):num_clus),
(sample(seq_len(num_clus / 2)))
)]
etaw <- rnorm(num_clus * clus_size, sd = psiw[clus_id])
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)Model Formulation
The model can be viewed as a conditional normal model. Specifically, if we condition on the cluster-specific random effects, the observations within a cluster follow a multivariate normal distribution.
Let \(\boldsymbol{\mathbf{\eta}}^b_j\) be the between-cluster latent factor and \(\log(\psi^w_j)\) be the log of the within-cluster factor variance. We also allow the factor loadings \(\boldsymbol{\mathbf{\Lambda}}_j\) to be random. Conditioned on these random effects, the within-cluster moments are:
\[ \begin{aligned} \bar{\boldsymbol{\mathbf{y}}}_j \mid \eta^b_j, \log(\psi^w_j), \boldsymbol{\mathbf{\Lambda}}_j & \sim N(\boldsymbol{\mathbf{\nu }}+ \boldsymbol{\mathbf{\Lambda}}_j \eta^b_j, \boldsymbol{\mathbf{\Sigma}}_{Wj} + \boldsymbol{\mathbf{\Theta}}^b) \\ \boldsymbol{\mathbf{y}}_{ij} - \bar{\boldsymbol{\mathbf{y}}}_j \mid \eta^b_j, \log(\psi^w_j), \boldsymbol{\mathbf{\Lambda}}_j & \sim N(\boldsymbol{\mathbf{0}}, \boldsymbol{\mathbf{\Sigma}}_{Wj}) \quad \text{for } i = 1, \ldots, n_j - 1 \end{aligned} \]
where the within-cluster covariance matrix is:
\[ \boldsymbol{\mathbf{\Sigma}}_{Wj} = \boldsymbol{\mathbf{\Lambda}}_j \psi^w_j \boldsymbol{\mathbf{\Lambda}}_j^\intercal + \boldsymbol{\mathbf{\Theta}}^w \]
Note that \(\psi^w_j = \exp(\log(\psi^w_j))\).
The between-subject random effects include the latent factor part \(\boldsymbol{\mathbf{b}}^\eta_j = (\eta^b_j, \log[\psi^w_j])\) and the item parameter part. For the items, we have random intercepts (which are standard in multilevel models) and random loadings. Let’s denote the vector of random loadings as \(\boldsymbol{\mathbf{b}}^\lambda_j = \mathop{\mathrm{vec}}(\boldsymbol{\mathbf{\Lambda}}_j)\) and random intercepts as \(\boldsymbol{\mathbf{\varepsilon}}^{b}_j\).
So assume we have four items and one factor, the full vector of random effects has a dimension of \(2 + 4 + 4 = 10\). A full \(10 \times 10\) covariance matrix is difficult to estimate. However, in factor analysis we usually assume:
- The latent factor components (\(\boldsymbol{\mathbf{b}}^\eta_j\)) are independent of the item-specific components (\(\boldsymbol{\mathbf{b}}^\lambda_j\) and \(\boldsymbol{\mathbf{\varepsilon}}^{b}_j\)).
- Local independence: The item-specific components for item \(k\) are independent of those for item \(k'\) for \(k \neq k'\).
Therefore, we only need to estimate: * The \(2 \times 2\) covariance matrix for the latent factor components (\(\eta^b_j\) and \(\log[\psi^w_j]\)). * The covariances between the random loading and random intercept for the same item (if both are random).
Using Mplus
Mplus allows for “random item parameters” using the | symbol in the %WITHIN% block. Here, lam1-lam4 | fw BY y1-y4 defines the random loadings. We then model these random variables in the %BETWEEN% block.
Code
# Save data for Mplus
write.table(dat, file = "mcfa_rl.dat", row.names = FALSE, col.names = FALSE)
writeLines(
"
TITLE: Bayesian 1-factor model with random variance and loadings;
DATA: FILE = mcfa_rl.dat;
VARIABLE:
NAMES = y1-y4 id;
USEVAR = y1-y4;
CLUSTER = id;
ANALYSIS:
TYPE = TWOLEVEL RANDOM;
ESTIMATOR = BAYES; ! This is default
BITERATION = 500000 (10000);
BCONVERGENCE = .005; ! ESS of ~ 100
CHAINS = 3;
PROCESS = 3;
MODEL: %WITHIN%
lam1-lam4 | fw BY y1-y4;
logpsiw | fw;
%BETWEEN%
fb BY y1-y4* (l1-l4);
logpsiw with fb; ! psi_12
fb; logpsiw; ! psi_1 and psi_2
logpsiw fb WITH lam1-lam4@0;
[logpsiw@0]; ! for setting the scale
[lam1-lam4] (l1-l4);
lam1-lam4*;
y1-y4;
lam1-lam4 PWITH y1-y4;
OUTPUT: TECH1 TECH8;
SAVEDATA:
BPARAMETERS = mcfa_random_load_draws.dat;
",
con = "mcfa_random_load.inp"
)Code
mplus_draws <- readModels("mcfa_random_load.out")$bparameters[[2]] # discard burn-in
# Drop chain and iteration columns
mplus_draws[[1]] <- mplus_draws[[1]][, -c(1, 2)]
mplus_draws[[2]] <- mplus_draws[[2]][, -c(1, 2)]
mplus_draws[[3]] <- mplus_draws[[3]][, -c(1, 2)]
mplus_draws <- posterior::as_draws_list(mplus_draws)
mplus_draws |>
posterior::summarise_draws() |>
knitr::kable(digits = 2)| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| Parameter.1_%WITHIN%:.Y1 | 0.53 | 0.53 | 0.01 | 0.01 | 0.51 | 0.55 | 1.00 | 7865.59 | 27926.49 |
| Parameter.2_%WITHIN%:.Y2 | 0.50 | 0.50 | 0.01 | 0.01 | 0.47 | 0.52 | 1.00 | 6027.18 | 21701.41 |
| Parameter.3_%WITHIN%:.Y3 | 0.40 | 0.40 | 0.01 | 0.01 | 0.39 | 0.42 | 1.00 | 9403.61 | 28093.05 |
| Parameter.4_%WITHIN%:.Y4 | 0.38 | 0.38 | 0.01 | 0.01 | 0.36 | 0.40 | 1.00 | 13012.29 | 30204.25 |
| Parameter.5_%BETWEEN%:.MEAN.LAM1(equality/label) | 0.63 | 0.62 | 0.02 | 0.02 | 0.60 | 0.66 | 1.02 | 112.20 | 246.34 |
| Parameter.6_%BETWEEN%:.MEAN.LAM2(equality/label) | 0.81 | 0.81 | 0.02 | 0.02 | 0.78 | 0.85 | 1.02 | 105.27 | 233.98 |
| Parameter.7_%BETWEEN%:.MEAN.LAM3(equality/label) | 0.62 | 0.62 | 0.02 | 0.02 | 0.59 | 0.66 | 1.02 | 111.53 | 299.16 |
| Parameter.8_%BETWEEN%:.MEAN.LAM4(equality/label) | 0.72 | 0.72 | 0.02 | 0.02 | 0.68 | 0.76 | 1.01 | 168.93 | 498.67 |
| Parameter.9_%BETWEEN%:.MEAN.Y1 | 0.02 | 0.02 | 0.02 | 0.02 | -0.01 | 0.05 | 1.02 | 144.79 | 497.37 |
| Parameter.10_%BETWEEN%:.MEAN.Y2 | 0.31 | 0.31 | 0.02 | 0.02 | 0.27 | 0.34 | 1.01 | 158.31 | 441.26 |
| Parameter.11_%BETWEEN%:.MEAN.Y3 | -0.49 | -0.49 | 0.02 | 0.02 | -0.53 | -0.45 | 1.01 | 322.34 | 1294.41 |
| Parameter.12_%BETWEEN%:.MEAN.Y4 | 0.20 | 0.20 | 0.02 | 0.02 | 0.17 | 0.24 | 1.01 | 233.03 | 785.25 |
| Parameter.13_%BETWEEN%:.FB | 0.14 | 0.14 | 0.02 | 0.02 | 0.11 | 0.18 | 1.00 | 458.60 | 1891.12 |
| Parameter.14_%BETWEEN%:.LAM1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 | 309.96 | 727.62 |
| Parameter.15_%BETWEEN%:.LAM2 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.01 | 266.85 | 846.06 |
| Parameter.16_%BETWEEN%:.LAM3 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.01 | 425.55 | 1339.94 |
| Parameter.17_%BETWEEN%:.LAM4 | 0.06 | 0.06 | 0.01 | 0.01 | 0.05 | 0.08 | 1.00 | 635.49 | 2543.96 |
| Parameter.18_%BETWEEN%:.LOGPSIW.WITH.FB | -0.26 | -0.26 | 0.03 | 0.03 | -0.32 | -0.21 | 1.00 | 889.12 | 3440.92 |
| Parameter.19_%BETWEEN%:.LOGPSIW | 0.89 | 0.89 | 0.10 | 0.09 | 0.74 | 1.06 | 1.01 | 732.80 | 2335.54 |
| Parameter.20_%BETWEEN%:.Y1.WITH.LAM1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.02 | 207.84 | 456.11 |
| Parameter.21_%BETWEEN%:.Y1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.02 | 141.58 | 222.45 |
| Parameter.22_%BETWEEN%:.Y2.WITH.LAM2 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.02 | 214.27 | 656.68 |
| Parameter.23_%BETWEEN%:.Y2 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00 | 0.01 | 1.02 | 174.99 | 206.09 |
| Parameter.24_%BETWEEN%:.Y3.WITH.LAM3 | 0.00 | 0.00 | 0.00 | 0.00 | -0.01 | 0.01 | 1.00 | 704.85 | 1721.15 |
| Parameter.25_%BETWEEN%:.Y3 | 0.11 | 0.11 | 0.01 | 0.01 | 0.09 | 0.12 | 1.00 | 16930.53 | 31074.11 |
| Parameter.26_%BETWEEN%:.Y4.WITH.LAM4 | -0.01 | -0.01 | 0.01 | 0.01 | -0.02 | 0.00 | 1.00 | 4120.75 | 11428.62 |
| Parameter.27_%BETWEEN%:.Y4 | 0.06 | 0.06 | 0.01 | 0.01 | 0.05 | 0.07 | 1.00 | 7310.45 | 13462.19 |
Using STAN
We will use a similar strategy as in Part II to efficiently sample from this model. The key is to recognize that constructing the full covariance matrix for all clusters in a straightforward way can be computationally expensive.
Instead, we can exploit the conditional structure. If we condition on the random loadings \(\boldsymbol{\mathbf{\Lambda}}_j\) and the random variance \(\psi^w_j\), the problem reduces to a more standard normal model which allows for sufficient statistics to be used. The Stan code below implements this by integrating the appropriate likelihoods.
Code
//
// This Stan program defines a one-factor two-level CFA model with a
// marginal likelihood approach similar to one described in
// https://ecmerkle.github.io/blavaan/articles/multilevel.html.
//
// Define function for computing multi-normal log density using sufficient
// statistics
functions {
real multi_normal_suff_lpdf(matrix S,
matrix Sigma,
int n) {
return -0.5 * (n * (log_determinant_spd(Sigma) + cols(S) * log(2 * pi())) +
trace(mdivide_left_spd(Sigma, S)));
}
}
// The input data are (a) sw, cluster-specific within-cluster cross-product
// matrix, (b) ybar, vectors of cluster means, and (c) n, cluster sample
// size. It assumes complete data.
data {
int<lower=1> p; // number of items
int<lower=1> J; // number of clusters
array[J] int<lower=1> n; // cluster size
array[J] matrix[p, p] sw; // within cross-product matrices
array[J] vector[p] ybar; // cluster means
}
// The parameters accepted by the model. Our model
// accepts four parameters:
// 'lambda_star' for unscaled loadings, 'nu' for intercepts,
// 'thetaw' and 'thetab' for uniqueness at the within and the between levels,
// and 'psib' for latent factor SD
parameters {
vector[p] mu_lambda_star; // unscaled mean loading vector assuming cross-level
// invariance
vector<lower=0>[p] thetaw; // within-cluster uniqueness vector
// vector<lower=0>[p] thetab; // between-cluster uniqueness vector
array[p] cholesky_factor_corr[2] Lthetab;
matrix<lower=0>[p, 2] sd_thetab; // SD of (dlambda, eps^b) for p items
cholesky_factor_corr[2] Lphi; // cholesky of correlation matrix of (logvar, mean)
real<lower=0> sd_logpsiw; // SD of log(SD) of etaw
real<lower=0> psib; // between-cluster latent standard deviation
row_vector[J] z_logpsiw; // scaled value of log(SD) of etaw
vector[p] nu; // item means/intercepts
matrix[p, J] z_dlambda; // scaled value of deviation from mean loading
}
// The model to be estimated. We compute the log-likelihood as the sum of
// two components: cluster-mean centered variables are multivariate normal
// with covariance matrix sigmaw (and degree of freedom adjustment),
// and cluster-mean variables are multivariate normal with covariance
// matrix sigmab + sigmaw / nj
model {
matrix[p, J] lambda = rep_matrix(mu_lambda_star, J) +
diag_pre_multiply(sd_thetab[, 1], z_dlambda);
row_vector[J] psiw = exp(sd_logpsiw .* z_logpsiw);
// Conditional value of thetab given dlambda;
matrix[p, 2] cond_thetab;
for (i in 1:p) {
Lthetab[i] ~ lkj_corr_cholesky(2);
cond_thetab[i, ] = Lthetab[i][2, ];
}
// Conditional mean given log(SD) of etaw and dlambda
matrix[p, J] mu = diag_post_multiply(lambda,
(psib * Lphi[2, 1]) .* z_logpsiw) +
diag_pre_multiply(sd_thetab[, 2] .* cond_thetab[, 1],
z_dlambda);
mu_lambda_star ~ normal(0, 10);
thetaw ~ gamma(1, .5);
// thetab ~ gamma(1, .5);
psib ~ gamma(1, .5);
nu ~ normal(0, 32);
z_logpsiw ~ std_normal();
to_vector(z_dlambda) ~ std_normal();
sd_logpsiw ~ gamma(1, 0.5);
Lphi ~ lkj_corr_cholesky(2);
for (j in 1:J) {
// The following becomes cluster-specific
matrix[p, p] LLt = lambda[, j] * lambda[, j]';
matrix[p, p] sigmab = add_diag(square(psib * Lphi[2, 2]) .* LLt,
square(sd_thetab[, 2] .* cond_thetab[, 2]));
matrix[p, p] sigmawj = add_diag(square(psiw[j]) .* LLt, square(thetaw));
target += multi_normal_suff_lpdf(sw[j] | sigmawj, n[j] - 1);
ybar[j] ~ multi_normal(nu + mu[, j], sigmab + sigmawj / n[j]);
}
}
// Rotating the factor loadings
// (see https://discourse.mc-stan.org/t/latent-factor-loadings/1483)
generated quantities {
vector[p] mu_lambda; // final loading vector
real cor_logpsiw_etab;
{
int sign_l1 = mu_lambda_star[1] > 0 ? 1 : -1;
mu_lambda = sign_l1 * mu_lambda_star;
cor_logpsiw_etab = sign_l1 * Lphi[2, 1];
}
}Code
# Prepare data for STAN
yc <- sweep(dat[1:4], MARGIN = 2, STATS = colMeans(dat[1:4])) # centered data
nj <- table(dat$id)
# Sufficient statistics
sswj <- tapply(yc, INDEX = dat$id, FUN = \(x) tcrossprod(t(x) - colMeans(x)))
ybarj <- tapply(dat[1:4], INDEX = dat$id, FUN = \(x) colMeans(x))Code
# Run STAN
mcfa_randload_fit <- mcfa_randload$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
)Running MCMC with 3 parallel chains...
Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup) Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup) Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup) Chain 3 Iteration: 500 / 1500 [ 33%] (Warmup)
Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
Chain 1 Iteration: 500 / 1500 [ 33%] (Warmup)
Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling) Chain 2 Iteration: 500 / 1500 [ 33%] (Warmup)
Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
Chain 3 Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain 1 Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain 2 Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
Chain 3 finished in 313.8 seconds.
Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
Chain 1 finished in 327.9 seconds.
Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
Chain 2 finished in 336.2 seconds.
All 3 chains finished successfully.
Mean chain execution time: 326.0 seconds.
Total execution time: 347.5 seconds.Code
mcfa_randload_fit$summary(
c(
"mu_lambda",
"thetaw",
"sd_thetab",
"sd_logpsiw",
"psib",
"cor_logpsiw_etab"
)
) |>
knitr::kable(digits = 2)| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| mu_lambda[1] | 0.62 | 0.62 | 0.02 | 0.02 | 0.59 | 0.65 | 1.00 | 722.18 | 1138.66 |
| mu_lambda[2] | 0.80 | 0.80 | 0.02 | 0.02 | 0.77 | 0.84 | 1.00 | 655.19 | 1059.48 |
| mu_lambda[3] | 0.62 | 0.62 | 0.02 | 0.02 | 0.59 | 0.65 | 1.00 | 699.90 | 1243.34 |
| mu_lambda[4] | 0.71 | 0.71 | 0.02 | 0.02 | 0.67 | 0.74 | 1.00 | 687.53 | 1048.89 |
| thetaw[1] | 0.73 | 0.73 | 0.01 | 0.01 | 0.71 | 0.74 | 1.00 | 3629.56 | 2275.06 |
| thetaw[2] | 0.71 | 0.71 | 0.01 | 0.01 | 0.69 | 0.73 | 1.00 | 3804.96 | 1666.89 |
| thetaw[3] | 0.64 | 0.64 | 0.01 | 0.01 | 0.62 | 0.65 | 1.00 | 3833.02 | 2167.82 |
| thetaw[4] | 0.62 | 0.62 | 0.01 | 0.01 | 0.60 | 0.63 | 1.01 | 2952.01 | 2159.50 |
| sd_thetab[1,1] | 0.05 | 0.05 | 0.02 | 0.02 | 0.01 | 0.08 | 1.01 | 362.74 | 846.65 |
| sd_thetab[2,1] | 0.04 | 0.04 | 0.02 | 0.03 | 0.00 | 0.08 | 1.01 | 384.40 | 637.47 |
| sd_thetab[3,1] | 0.02 | 0.02 | 0.02 | 0.01 | 0.00 | 0.05 | 1.00 | 461.68 | 961.37 |
| sd_thetab[4,1] | 0.24 | 0.24 | 0.02 | 0.02 | 0.21 | 0.27 | 1.00 | 858.77 | 1654.67 |
| sd_thetab[1,2] | 0.03 | 0.03 | 0.02 | 0.02 | 0.00 | 0.07 | 1.00 | 1502.14 | 1313.26 |
| sd_thetab[2,2] | 0.05 | 0.05 | 0.03 | 0.03 | 0.00 | 0.10 | 1.00 | 1162.27 | 1513.45 |
| sd_thetab[3,2] | 0.32 | 0.32 | 0.02 | 0.02 | 0.30 | 0.35 | 1.00 | 4011.86 | 2062.37 |
| sd_thetab[4,2] | 0.23 | 0.23 | 0.02 | 0.02 | 0.21 | 0.26 | 1.00 | 4040.65 | 2309.13 |
| sd_logpsiw | 0.47 | 0.47 | 0.03 | 0.03 | 0.43 | 0.51 | 1.00 | 1235.37 | 1782.72 |
| psib | 0.39 | 0.39 | 0.03 | 0.03 | 0.34 | 0.43 | 1.00 | 1585.03 | 1892.98 |
| cor_logpsiw_etab | -0.72 | -0.72 | 0.06 | 0.06 | -0.81 | -0.61 | 1.00 | 1488.01 | 1947.50 |
Code
data.frame(
model = c("Mplus (Gibbs)", "STAN"),
ess = c(
posterior::summarise_draws(mplus_draws, "ess_bulk")[[2]][21],
mcfa_randload_fit$summary(
variables = "sd_thetab[4,1]",
posterior::ess_bulk
)[[2]]
),
time = c(
mplus_sec,
mcfa_randload_fit$time()$total
)
) |>
knitr::kable(
digits = 2,
col.names = c(
"Model",
"ESS (Bulk) for the variance/SD of $b^\\lambda_4$",
"Time (seconds)"
)
)| Model | ESS (Bulk) for the variance/SD of \(b^\lambda_4\) | Time (seconds) |
|---|---|---|
| Mplus (Gibbs) | 141.58 | 1152.0 |
| STAN | 858.77 | 347.5 |
The sampling is much more efficient for STAN than for Mplus in this model, as the ESS is larger with the STAN approach for most parameters.
NoteSession Info
Code
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.2 (2025-10-31)
os Ubuntu 24.04.4 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate C.UTF-8
ctype C.UTF-8
tz Asia/Macau
date 2026-02-13
pandoc 3.6.3 @ /home/marklai/.positron-server/bin/4b472e9bf40ca88f886fe9c7fcfbb9ee702fccd0/quarto/bin/tools/x86_64/ (via rmarkdown)
quarto 1.8.27 @ /home/marklai/.positron-server/bin/4b472e9bf40ca88f886fe9c7fcfbb9ee702fccd0/quarto/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] RSPM (R 4.5.0)
backports 1.5.0 2024-05-23 [1] RSPM (R 4.5.0)
boot 1.3-32 2025-08-29 [4] CRAN (R 4.5.2)
checkmate 2.3.4 2026-02-03 [1] RSPM (R 4.5.0)
cli 3.6.5 2025-04-23 [1] RSPM (R 4.5.0)
cmdstanr * 0.9.0 2025-03-30 [1] https://stan-dev.r-universe.dev (R 4.5.2)
coda 0.19-4.1 2024-01-31 [1] RSPM (R 4.5.0)
data.table 1.18.2.1 2026-01-27 [1] RSPM (R 4.5.0)
digest 0.6.39 2025-11-19 [1] RSPM (R 4.5.0)
distributional 0.6.0 2026-01-14 [1] RSPM (R 4.5.0)
dplyr 1.2.0 2026-02-03 [1] RSPM (R 4.5.0)
evaluate 1.0.5 2025-08-27 [1] RSPM (R 4.5.0)
farver 2.1.2 2024-05-13 [1] RSPM (R 4.5.0)
fastDummies 1.7.5 2025-01-20 [1] RSPM (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] RSPM (R 4.5.0)
generics 0.1.4 2025-05-09 [1] RSPM (R 4.5.0)
ggplot2 4.0.2 2026-02-03 [1] RSPM (R 4.5.0)
glue 1.8.0 2024-09-30 [1] RSPM (R 4.5.0)
gsubfn 0.7 2018-03-16 [1] RSPM (R 4.5.0)
gtable 0.3.6 2024-10-25 [1] RSPM (R 4.5.0)
htmltools 0.5.9 2025-12-04 [1] RSPM (R 4.5.0)
htmlwidgets 1.6.4 2023-12-06 [1] RSPM (R 4.5.0)
httr 1.4.7 2023-08-15 [1] RSPM (R 4.5.0)
jsonlite 2.0.0 2025-03-27 [1] RSPM (R 4.5.0)
knitr 1.51 2025-12-20 [1] RSPM (R 4.5.0)
lattice 0.22-7 2025-04-02 [4] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] RSPM (R 4.5.0)
magrittr 2.0.4 2025-09-12 [1] RSPM (R 4.5.0)
matrixStats 1.5.0 2025-01-07 [1] RSPM (R 4.5.0)
MplusAutomation * 1.2 2025-09-02 [1] RSPM (R 4.5.0)
otel 0.2.0 2025-08-29 [1] RSPM (R 4.5.0)
pander 0.6.6 2025-03-01 [1] RSPM (R 4.5.0)
pillar 1.11.1 2025-09-17 [1] RSPM (R 4.5.0)
pkgconfig 2.0.3 2019-09-22 [1] RSPM (R 4.5.0)
plyr 1.8.9 2023-10-02 [1] RSPM (R 4.5.0)
posterior * 1.6.1.9000 2026-02-01 [1] https://stan-dev.r-universe.dev (R 4.5.2)
processx 3.8.6 2025-02-21 [1] RSPM (R 4.5.0)
proto 1.0.0 2016-10-29 [1] RSPM (R 4.5.0)
ps 1.9.1 2025-04-12 [1] RSPM (R 4.5.0)
R6 2.6.1 2025-02-15 [1] RSPM (R 4.5.0)
RColorBrewer 1.1-3 2022-04-03 [1] RSPM (R 4.5.0)
Rcpp 1.1.1 2026-01-10 [1] RSPM (R 4.5.0)
rlang 1.1.7 2026-01-09 [1] RSPM (R 4.5.0)
rmarkdown 2.30 2025-09-28 [1] RSPM (R 4.5.0)
S7 0.2.1 2025-11-14 [1] RSPM (R 4.5.0)
scales 1.4.0 2025-04-24 [1] RSPM (R 4.5.0)
sessioninfo 1.2.3 2025-02-05 [1] RSPM (R 4.5.0)
tensorA 0.36.2.1 2023-12-13 [1] RSPM (R 4.5.0)
texreg 1.39.5 2025-12-22 [1] RSPM (R 4.5.0)
tibble 3.3.1 2026-01-11 [1] RSPM (R 4.5.0)
tidyselect 1.2.1 2024-03-11 [1] RSPM (R 4.5.0)
vctrs 0.7.1 2026-01-23 [1] RSPM (R 4.5.0)
withr 3.0.2 2024-10-28 [1] RSPM (R 4.5.0)
xfun 0.56 2026-01-18 [1] RSPM (R 4.5.0)
xtable 1.8-4 2019-04-21 [1] RSPM (R 4.5.0)
yaml 2.3.12 2025-12-10 [1] RSPM (R 4.5.0)
[1] /home/marklai/R/x86_64-pc-linux-gnu-library/4.5
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────