EFA example based off Rick Farouni's example on Holzinger-Swineford data

0_efa_example.R

library(rstan)
library(lavaan)

cmdstanr::set_cmdstan_path("~/cmdstan/")

# https://rfarouni.github.io/assets/projects/BayesianFactorAnalysis/BayesianFactorAnalysis.html
# L is N_items BY N_factor loading matrix, constrained lower-triangular
# Cov = LL' + diag(res_vars), factors constrained orthogonal

# Chose items 3, 6 and 8 as markers of factors 1, 2 and 3
# I expect them to have positive loadings + high local independence
dat_list <- list(
  input_data = HolzingerSwineford1939[, paste0("x", 1:9)],  # raw data
  markers = c(3, 6, 8),  # items to identify each factor
  Nf = 3  # number of factors
)
dat_list$S <- cov(dat_list$input_data)
dat_list$Np <- nrow(dat_list$input_data)  # number persons
dat_list$Ni <- ncol(dat_list$input_data)  # number items
dat_list

# Raw data approach
(efa.stan <- cmdstanr::cmdstan_model(file.path("efa_farouni.stan")))
# Working off sample covariance matrix, much faster
(efa_cov.stan <- cmdstanr::cmdstan_model(file.path("efa_cov_farouni.stan")))
efa_0 <- efa_cov.stan$sample(
  data = dat_list, seed = 12345, iter_warmup = 1e3, iter_sampling = 1e3,
  chains = 3, parallel_chains = 4, refresh = 200)
# Warnings on L should be ignored because of sign-switching,
# loading_matrix parameter is corrected varion of L
efa_0$cmdstan_diagnose()
efa_0 <- read_stan_csv(efa_0$output_files())
print(efa_0, c("intercepts", "L"), include = FALSE, digits = 4)
stan_hist(efa_0, c("intercepts", "L"), include = FALSE)
traceplot(efa_0, c("intercepts", "L"), include = FALSE)

# Kind of matches typical Holzinger-Swineford solution
matrix(colMeans(as.data.frame(efa_0, "loading_matrix")), dat_list$Ni)

efa_cov_farouni.stan

data {
  int<lower = 0> Np; // N_persons
  int<lower = 0> Ni; // N_items
  int<lower = 0> Nf; // N_factors
  cov_matrix[Ni] S; // input matrix
  int markers[Nf]; // markers for factors
}
transformed data {
  int Nl = 0;  // number of loadings
  int loading_pattern[Ni, Nf] = rep_array(1, Ni, Nf);

  for (i in 1:Nf) Nl += Ni - i + 1;

  for (j in 1:(Nf - 1)) {
    for (k in (j + 1):Nf) {
      loading_pattern[markers[j], k] = 0;
    }
  }
}
parameters {
  real<lower = 0> l_sd;  // parameter to constrain loadings except diagonal of loading mat
  vector[Nl] L;  // loadings
  vector<lower = 0>[Ni] item_res_sd;  // residual SDs
}
model {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0.0, Ni, Nf);

  {
    int pos = 0;
    for (j in 1:Nf) {
      for (i in 1:Ni) {
        if (loading_pattern[i, j] == 1) {
          pos += 1;
          loading_matrix[i, j] = L[pos] * l_sd;
          if (i == markers[j]) loading_matrix[i, j] /= l_sd;
        }
      }
    }
  }

  l_sd ~ std_normal();
  L ~ std_normal();
  item_res_sd ~ student_t(3, 0, 1);

  {
    matrix[Ni, Ni] Sigma;

    Sigma = multiply_lower_tri_self_transpose(loading_matrix) +
      diag_matrix(square(item_res_sd));
    S ~ wishart(Np - 1, Sigma / (Np - 1));
  }
}
generated quantities {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0.0, Ni, Nf);

  {
    int pos = 0;
    for (j in 1:Nf) {
      for (i in 1:Ni) {
        if (loading_pattern[i, j] == 1) {
          pos += 1;
          loading_matrix[i, j] = L[pos] * l_sd;
          if (i == markers[j]) loading_matrix[i, j] /= l_sd;
        }
      }
      if (loading_matrix[markers[j], j] < 0) {
        loading_matrix[, j] *= -1;
      }
    }
  }
}

efa_farouni.stan

data {
  int<lower = 0> Np; // N_persons
  int<lower = 0> Ni; // N_items
  int<lower = 0> Nf; // N_factors
  vector[Ni] input_data[Np]; // input matrix
  int markers[Nf]; // markers for factors
}
transformed data {
  int Nl = 0;  // number of loadings
  int loading_pattern[Ni, Nf] = rep_array(1, Ni, Nf);

  for (i in 1:Nf) Nl += Ni - i + 1;

  for (j in 1:(Nf - 1)) {
    for (k in (j + 1):Nf) {
      loading_pattern[markers[j], k] = 0;
    }
  }
}
parameters {
  real<lower = 0> l_sd;  // parameter to constrain loadings except diagonal of loading mat
  vector[Nl] L;  // loadings
  vector<lower = 0>[Ni] item_res_sd;  // residual SDs
  vector[Ni] intercepts;  // intercepts
}
model {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0.0, Ni, Nf);

  {
    int pos = 0;
    for (j in 1:Nf) {
      for (i in 1:Ni) {
        if (loading_pattern[i, j] == 1) {
          pos += 1;
          loading_matrix[i, j] = L[pos] * l_sd;
          if (i == markers[j]) loading_matrix[i, j] /= l_sd;
        }
      }
    }
  }

  l_sd ~ std_normal();
  L ~ std_normal();
  item_res_sd ~ student_t(3, 0, 1);
  intercepts ~ normal(0, 10);

  {
    matrix[Ni, Ni] Sigma;

    Sigma = multiply_lower_tri_self_transpose(loading_matrix) +
      diag_matrix(square(item_res_sd));
    input_data ~ multi_normal(intercepts, Sigma);
  }
}
generated quantities {
  matrix[Ni, Nf] loading_matrix = rep_matrix(0.0, Ni, Nf);

  {
    int pos = 0;
    for (j in 1:Nf) {
      for (i in 1:Ni) {
        if (loading_pattern[i, j] == 1) {
          pos += 1;
          loading_matrix[i, j] = L[pos] * l_sd;
          if (i == markers[j]) loading_matrix[i, j] /= l_sd;
        }
      }
      if (loading_matrix[markers[j], j] < 0) {
        loading_matrix[, j] *= -1;
      }
    }
  }
}

Based off: https://rfarouni.github.io/assets/projects/BayesianFactorAnalysis/BayesianFactorAnalysis.html