Holzinger Swineford Bayesian CFA example

0_hs_gaus.R

library(lavaan)

summary(m0 <- cfa(
  "F1 =~ x1 + x2 + x3 + x9\n F2 =~ x4 + x5 + x6\n F3 =~ x7 + x8 + x9\n x3 ~~ x5\n x2 ~~ x7\n x4 ~~ x7",
  HolzingerSwineford1939, std.lv = TRUE, meanstructure = TRUE))

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# Cov mat decomp
# saveRDS(stan_model(stanc_ret = stanc(file = "stan_scripts/cfa_decomp.stan")),
#         "stan_scripts/cfa_decomp.rds")
cfa.decomp.stan <- readRDS("stan_scripts/cfa_decomp.rds")
# Cov mat decomp w/ residual covariances
# saveRDS(stan_model(stanc_ret = stanc(file = "stan_scripts/cfa_decomp_pe.stan")),
#         "stan_scripts/cfa_decomp_pe.rds")
cfa.decomp.pe.stan <- readRDS("stan_scripts/cfa_decomp_pe.rds")

X <- HolzingerSwineford1939[, paste0("x", 1:9)]

(loading_pat <- matrix(c(
  rep(1, 3), rep(0, 5), 1, rep(0, 3), rep(1, 3), rep(0, 9), rep(1, 3)), 9, 3))
dat <- list(
  Np = nrow(X), # number respondents
  Ni = ncol(X), # number items
  Nf = ncol(loading_pat), # number factors
  Nl = sum(loading_pat), # number loadings
  input_data = X, # input data matrix
  loading_pattern = loading_pat, # Loading patterns
  g_alpha = 1, g_beta = 1, # Inverse gamma priors for residual variances
  shape_phi_c = 2, # LKJ prior for interfactor correlation matrix
  # Column matrix with items containing items with loadings constrained to be positive
  # Useful for identifying the factors
  markers = matrix(c(1, 4, 7)),
  # Remaining items are only relevant if you want to model residual covariances
  v = 3, # number residual covariances
  # matrix of residual covariances
  # each row contains a pair and the numbers are the item positions in the input data matrix
  error_mat = matrix(c(3, 5, 2, 7, 4, 7), 3, byrow = T),
  beta_par = 2 # Symmetric beta prior for residual correlations
)
# Initialize all loadings to 1
initf <- function() list(loadings_mark = rep(1, dat$Nf), loadings = rep(1, dat$Nl - dat$Nf))

fit.decomp <- sampling(cfa.decomp.stan, data = dat, init = initf)
fit.decomp.pe <- sampling(cfa.decomp.pe.stan, data = dat, init = initf)

print(fit.decomp, pars = c(
  "loadings_mark", "loadings", "item_res_vars", "intercepts",
  "phi_mat[1, 2]", "phi_mat[1, 3]", "phi_mat[2, 3]", "sigma_loadings"),
  probs = c(.025, .975), digits_summary = 3)
stan_hist(fit.decomp, pars = c(
  "loadings_mark", "loadings", "item_res_vars", "intercepts",
  "phi_mat[1, 2]", "phi_mat[1, 3]", "phi_mat[2, 3]", "sigma_loadings"))
print(fit.decomp.pe, pars = c(
  "loadings_mark", "loadings", "res_cov", "item_res_vars", "intercepts",
  "phi_mat[1, 2]", "phi_mat[1, 3]", "phi_mat[2, 3]", "sigma_loadings"),
  probs = c(.025, .975), digits_summary = 3)
stan_hist(fit.decomp.pe, pars = c(
  "loadings_mark", "loadings", "res_cov", "item_res_vars", "intercepts",
  "phi_mat[1, 2]", "phi_mat[1, 3]", "phi_mat[2, 3]", "sigma_loadings"))

library(loo)

(hs <- loo(fit.decomp))
(hs.pe <- loo(fit.decomp.pe))
compare(hs, hs.pe)
loo_model_weights(list(hs, hs.pe))

cfa_decomp.stan

// No residual covariances
data {
  real g_alpha; // prior for inv_gamma
  real g_beta; // prior for inv_gamma
  int<lower = 1> shape_phi_c; // lkj prior shape for phi
  int<lower = 0> Np; // N_persons
  int<lower = 0> Ni; // N_items
  int<lower = 0> Nf; // N_factors
  int<lower = 0> Nl; // N_non-zero loadings
  matrix[Np, Ni] input_data; // input matrix
  matrix[Ni, Nf] loading_pattern;
  int markers[Nf, 1]; // markers
}
transformed data {
  vector[Ni] dat_mv[Np];

  for (i in 1:Np) {
    dat_mv[i] = input_data[i, ]';
  }
}
parameters {
  vector<lower = 0>[Nf] loadings_mark; // loadings
  vector[Nl - Nf] loadings; // loadings
  real<lower = 0> sigma_loadings; // sd of loadings, hyperparm
  vector<lower = 0>[Ni] item_res_vars; // item residual vars heteroskedastic
  cholesky_factor_corr[Nf] phi_mat_chol;
  vector[Ni] intercepts;
}
transformed parameters {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0, Ni, Nf);
  corr_matrix[Nf] phi_mat;
  matrix[Ni, Ni] lamb_phi_lamb;
  matrix[Ni, Ni] delta_mat;
  matrix[Ni, Ni] Sigma;

  {
    int pos = 0;
    for (i in 1:Ni) {
      for (j in 1:Nf) {
        if (loading_pattern[i, j] != 0) {
          if (i == markers[j, 1]) {
            loading_matrix[i, j] = loadings_mark[j];
          } else {
            pos = pos + 1;
            loading_matrix[i, j] = loadings[pos];
          }
        }
      }
    }
  }

  phi_mat = multiply_lower_tri_self_transpose(phi_mat_chol);

  lamb_phi_lamb = quad_form_sym(phi_mat, loading_matrix');
  delta_mat = diag_matrix(item_res_vars);
  Sigma = lamb_phi_lamb + delta_mat;
}
model {
  loadings ~ normal(0, sigma_loadings);
  loadings_mark ~ normal(0, sigma_loadings);
  sigma_loadings ~ student_t(3, 0, 1);
  item_res_vars ~ inv_gamma(g_alpha, g_beta);
  phi_mat_chol ~ lkj_corr_cholesky(shape_phi_c);
  intercepts ~ normal(0, 10);

  dat_mv ~ multi_normal(intercepts, Sigma);
}
generated quantities {
  vector[Np] log_lik;

  for (i in 1:Np) {
    log_lik[i] = multi_normal_lpdf(dat_mv[i] | intercepts, Sigma);
  }
}

cfa_decomp_pe.stan

// Residual covariances via parameter expansion
functions {
  int sign(real x) {
    if (x > 0)
      return 1;
    else
      return -1;
  }
}
data {
  real g_alpha; // prior for inv_gamma
  real g_beta; // prior for inv_gamma
  real beta_par; // prior for correlations
  int<lower = 1> shape_phi_c; // lkj prior shape for phi
  int<lower = 0> Np; // N_persons
  int<lower = 0> Ni; // N_items
  int<lower = 0> Nf; // N_factors
  int<lower = 0> Nl; // N_non-zero loadings
  int<lower = 1> v; // N_correlated errors
  int error_mat[v, 2]; // cor error matrix
  matrix[Np, Ni] input_data; // input matrix
  matrix[Ni, Nf] loading_pattern;
  int markers[Nf, 1]; // markers
}
transformed data {
  vector[Ni] dat_mv[Np];

  for (i in 1:Np) {
    dat_mv[i] = input_data[i, ]';
  }
}
parameters {
  vector<lower = 0>[Nf] loadings_mark; // loadings
  vector[Nl - Nf] loadings; // loadings
  real<lower = 0> sigma_loadings; // sd of loadings, hyperparm
  vector<lower = 0>[Ni] item_res_vars; // item residual vars heteroskedastic
  cholesky_factor_corr[Nf] phi_mat_chol;
  vector[Ni] intercepts;
  vector<lower = 0, upper = 1>[v] cor_errs_01; // correlated errors on 01
}
transformed parameters {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0, Ni, Nf);
  matrix[Nf, Nf] phi_mat;
  matrix[Ni, Ni] lamb_phi_lamb;
  matrix[Ni, Ni] delta_mat_ast;
  matrix[Ni, Ni] delta_mat;
  matrix[Ni, Ni] Sigma;
  matrix[Ni, v] loading_par_exp;
  vector<lower = -1, upper = 1>[v] cor_errs; // correlated errors
  vector[v] res_cov; // residual covariances

  {
    int pos = 0;
    for (i in 1:Ni) {
      for (j in 1:Nf) {
        if (loading_pattern[i, j] != 0) {
          if (i == markers[j, 1]) {
            loading_matrix[i, j] = loadings_mark[j];
          } else {
            pos = pos + 1;
            loading_matrix[i, j] = loadings[pos];
          }
        }
      }
    }
  }

  phi_mat = multiply_lower_tri_self_transpose(phi_mat_chol);
  lamb_phi_lamb = quad_form_sym(phi_mat, loading_matrix');

  delta_mat = diag_matrix(item_res_vars);

  cor_errs = cor_errs_01 * 2 - 1;

  {
    int pos = 0;
    loading_par_exp = rep_matrix(0, Ni, v);
    for (i in 1:v) {
      res_cov[i] = cor_errs[i] * sqrt(prod(item_res_vars[error_mat[i, ]]));
      loading_par_exp[error_mat[i, 1], i] = sqrt(
        fabs(cor_errs[i]) * item_res_vars[error_mat[i, 1]]);
      loading_par_exp[error_mat[i, 2], i] = sign(cor_errs[i]) * sqrt(
        fabs(cor_errs[i]) * item_res_vars[error_mat[i, 2]]);
    }
  }

  {
    matrix[Ni, Ni] loading_par_exp_2;
    loading_par_exp_2 = tcrossprod(loading_par_exp);
    delta_mat_ast = diag_matrix(item_res_vars - diagonal(loading_par_exp_2));
    Sigma = lamb_phi_lamb + loading_par_exp_2 + delta_mat_ast;
  }

}
model {
  loadings ~ normal(0, sigma_loadings);
  loadings_mark ~ normal(0, sigma_loadings);
  sigma_loadings ~ student_t(3, 0, 1);
  item_res_vars ~ inv_gamma(g_alpha, g_beta);
  phi_mat_chol ~ lkj_corr_cholesky(shape_phi_c);
  intercepts ~ normal(0, 10);
  cor_errs_01 ~ beta(beta_par, beta_par);

  dat_mv ~ multi_normal(intercepts, Sigma);
}
generated quantities {
  vector[Np] log_lik;

  for (i in 1:Np) {
    log_lik[i] = multi_normal_lpdf(dat_mv[i] | intercepts, Sigma);
  }
}

Feel free to acknowledge, if you prefer to cite:

Uanhoro, J. O. (2019). Holzinger Swineford Bayesian CFA example. Retrieved ..., from https://gist.github.com/stonegold546/bc34a68e5cb338b05335c38273839bf0