# Multi Normal model

Hi all,

I am new to Stan, and I am trying to fit a CFA model with multivariate normal indicator variables. I am using Stan directly instead of blavaan because I am going to build on this model to include a hierarchical part to it. So the CFA part is just a start.

So these are the dimensions of my data.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

``````data {
int N; //number of persons
int ind; //number of indicator variables i.e. 9 with 3 variables per factor
matrix[N, ind] Y; // my input data
int Nf; //number of factors which will be 3
}

parameters {
matrix[2, Nf] ld; //only 2 loadings per factor since a reference variable will be defined for each factor
vector[ind] intercepts;
cholesky_factor_corr[Nf] Nf_corr;
vector < lower = 0 > [ind] sd_res; //residual sds
}

transformed parameters {
matrix[N, ind] mu_Y; //expected values of Y
matrix[N, Nf] LVs; //latent variables
vector[Nf] zeros;
cov_matrix[ind] Sigma;

zeros = rep_vector(0, Nf)
}

model {
//priors
for (k in 1:ind) {
intercepts[k] ~ normal(0, 10);
sd_res [k] ~ cauchy(0, 2.5);
}

for ( i in 1:2) {
for(j in 1:Nf) {
ld[i,j] ~ normal(0,10);
}
}

// distribution of data
to_vector(Y) ~ multi_normal(to_vector(mu_Y), Sigma);

// distribution of factors
to_vector(LVs) ~ multi_normal(zeros, Nf_corr);

}
``````

I keep getting errors based on the multi_normal_lpdf such as: Size of random variable (2709) and rows of covariance parameter (9) must match in size.
I have also left out some part of the code just to keep it simple.
I will be appreciate any help offered or if any one can point me in the direction of a STAN model in hierarchical SEM.

Thanks

It may help you to look at page 34 of the Stan User’s Guide for version 2.21. There is a model there that uses a multivariate normal, though for continuous predictors.

``````data {
int<lower=1> K;
int<lower=1> J;
int<lower=0> N;
vector[J] x[N];
vector[K] y[N];
}
parameters {
matrix[K, J] beta;
cov_matrix[K] Sigma;
}
model {
vector[K] mu[N];
for (n in 1:N)
mu[n] = beta*x[n];
y ~ multi_normal(mu, Sigma);
}
``````

In that case there are N observations of y and x. Each y is a vector of K elements and so y is declared as

``````vector[K] y[N];
``````

Each x is a vector of J predictors and x is declared as

``````vector[J] x[N];
``````

The coefficients beta that are applied to x to predict the mean y values is a K x J matrix

``````matrix[K,J] beta;
``````

The predicted mean values, mu, have the same kind of declaration as y

``````vector[K] mu[N];
``````

A for loop is used to generate the N vectors of mu and then the y values are predicted with

``````y ~ multi_normal(mu, Sigma);
``````

where Sigma is a K x K covariance matrix. That is vectorized and does not show explicitly that each y[n] which is a K vector is being predicted using a K vector of mu values. In your code, you have changed the Y and mu values into single vectors, which is not what multi-normal is expecting. Because your Y has 2709 values, it expects Sigma to be a 2709 x 2709 matrix. You need your Y to be an array of N vectors each with 9 values.

1 Like

Thank you very much, Fjcc!

I guess you got your problem solved. Would you mind posting the solution? I am curious to see what worked.

Yes, it runs for now without errors. I am still tweaking it though.

``````data {
int<lower = 1> Nf; // number of LVs
int<lower = 1> N; // number of persons
int <lower = 1> ind; //number of indicator variables
vector[ind] Y[N]; //data i.e. Y has a size of N, but each Y is a vector of 'ind' elements
}

parameters {
matrix[2,Nf] Nf_ld;
vector[ind] intercepts;
cholesky_factor_corr[Nf] Nf_corr;
vector <lower = 0> [ind] sd_res;
}

transformed parameters {

row_vector[Nf] ones;
matrix[3, Nf] ld_matrix;
cov_matrix[ind] Sigma;
cov_matrix[ind] theta_eps;
corr_matrix[Nf] omega;
matrix[ind,ind] lamba_phi_lampha;
vector[Nf] zeros;

zeros = rep_vector(0, Nf);

ones = rep_row_vector(1, Nf);
ld_matrix = append_row(ones, Nf_ld);

{
int pos = 0;
for (i in 1:ind){
for (j in 1:Nf) {

if(l_p[i,j] != 0) {
pos = pos + 1;
} else {
}
}
}
}

omega = multiply_lower_tri_self_transpose(Nf_corr);

theta_eps = diag_matrix(sd_res);

Sigma = lamba_phi_lampha + theta_eps;
}

model{
//vector[ind] mu_Y[N];
vector[Nf] LVs[N];

//distribution of factors
//LVs ~ multi_normal(zeros, Nf_corr);

//distribution of data
Y ~ multi_normal(intercepts, Sigma);

//priors
intercepts ~ normal(0,10);
sd_res ~ cauchy(0,3);

to_vector(Nf_ld) ~ normal(0,10);

Nf_corr ~ lkj_corr_cholesky(1);
}
``````
1 Like