Last active
August 16, 2019 13:57
-
-
Save stonegold546/bc34a68e5cb338b05335c38273839bf0 to your computer and use it in GitHub Desktop.
Holzinger Swineford Bayesian CFA example
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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); | |
} | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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