I started with some code from Mike Laurence and Michael de Witt that I modified to include an endogenous latent variable. I want to now make this multi-group by time so I can make the parameter factor loadings evolve over time. I am unsure how I should best incorporate multiple groups in this model. After that I think I should be able to add a dynamic component with a state-space model.
I got to the cholesky_factor_corr parameter declaration and got stumped at how to proceed vectorizing it for time/group. Any suggestions on how to perform this analysis by groups would be highly appreciated.
//confirmatory factor analysis code
data{
// n_subj: number of subjects
int n_subj ;
// n_y: number of outcomes
int n_y ;
// y: matrix of outcomes
matrix[n_subj,n_y] y ;
// n_fac: number of latent factors level one
int n_fac_one ;
// y_fac: list of which factor is associated with each outcome
int<lower=1,upper=n_fac_one> y_fac_one[n_y] ;
// n_fac: number of latent factors level two
int n_fac_two ;
// y_fac: list of which factor is associated with each outcome
int<lower=1,upper=n_fac_two> y_fac_two[n_fac_one] ;
}
transformed data{
// scaled_y: z-scaled outcomes
matrix[n_subj,n_y] scaled_y ;
for(i in 1:n_y){
scaled_y[,i] = (y[,i]-mean(y[,i])) / sd(y[,i]) ;
}
}
parameters{
// normal01: a helper variable for more efficient non-centered multivariate
// sampling of subj_facs
matrix[n_fac_one, n_subj] normal01;
matrix[n_fac_two, n_subj] normal02;
// fac_cor_helper: correlations (on cholesky factor scale) amongst latent
// factors
cholesky_factor_corr[n_fac_one] fac_cor_helper_one ;
cholesky_factor_corr[n_fac_two] fac_cor_helper_two ;
// scaled_y_means: means for each outcome
vector[n_y] scaled_y_means ;
// scaled_y_noise: measurement noise for each outcome
vector<lower=0>[n_y] scaled_y_noise_one ;
vector<lower=0>[n_y] scaled_y_noise_two ;
// betas: how each factor loads on to each outcome.
// Must be postive for identifiability.
vector<lower=0>[n_y] betas_one ;
vector<lower=0>[n_y] betas_two ;
}
transformed parameters{
// subj_facs: matrix of values for each factor associated with each subject
matrix[n_subj,n_fac_one] subj_facs_one ;
matrix[n_subj,n_fac_two] subj_facs_two ;
// the following calculation implies that rows of subj_facs are sampled from
// a multivariate normal distribution with mean of 0 and sd of 1 in all
// dimensions and a correlation matrix of fac_cor.
subj_facs_one = transpose(fac_cor_helper_one * normal01) ;
subj_facs_two = transpose(fac_cor_helper_two * normal02) ;
}
model{
// normal01 must have normal(0,1) prior for "non-centered" parameterization
// on the multivariate distribution of latent factors
to_vector(normal01) ~ normal(0,1) ;
to_vector(normal02) ~ normal(0,1) ;
// flat prior on correlations
fac_cor_helper_one ~ lkj_corr_cholesky(1) ;
fac_cor_helper_two ~ lkj_corr_cholesky(1) ;
// normal prior on y means
scaled_y_means ~ normal(0,1) ;
// weibull prior on measurement noise
scaled_y_noise_one ~ weibull(2,1) ;
scaled_y_noise_two ~ weibull(2,1) ;
// normal(0,1) priors on betas
betas_one ~ normal(0,1) ;
betas_two ~ normal(0,1) ;
// each outcome as normal, influenced by its associated latent factor
for(i in 1:n_y){
scaled_y[,i] ~ normal(
scaled_y_means[i] + subj_facs_one[,y_fac_one[i]] * betas_one[i]
, scaled_y_noise_one[i]
) ;
}
// each outcome as normal, influenced by its associated latent factor level two
for(i in 1:n_fac_one){
subj_facs_one[,i] ~ normal(
subj_facs_two[,y_fac_two[i]] * betas_two[i]
, scaled_y_noise_two[i]
) ;
}
}
generated quantities{
//cors: correlations amongst within-subject predictors' effects
corr_matrix[n_fac_one] fac_cor_one ;
corr_matrix[n_fac_two] fac_cor_two ;
// y_means: outcome means on the original scale of Y
vector[n_y] y_means ;
// y_noise: outcome noise on the original scale of Y
vector[n_y] y_noise ;
fac_cor_one = multiply_lower_tri_self_transpose(fac_cor_helper_one) ;
fac_cor_two = multiply_lower_tri_self_transpose(fac_cor_helper_two) ;
for(i in 1:n_y){
y_means[i] = scaled_y_means[i]*sd(y[,i]) + mean(y[,i]) ;
y_noise[i] = scaled_y_noise_one[i]*sd(y[,i]) ;
}
}
instead of
model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}
To include mathematical notation in your post put LaTeX syntax between two $ symbols, e.g.,
p(\theta | y) \propto p(\theta) p(y | \theta).
