I am fitting a simple, one-factor latent variable model in Stan, comparing the results to the output of `lavaan`

, but am running into problems estimating models with different identification constraints. In particular, setting the latent variable to be a unit normal and estimating all factor loadings and intercepts.

Consider J measured variables y_{ij} for instances i = 1 ... N and variables j = 1 ... J, each with their own intercept parameter \alpha_j and factor loading \beta_j on latent variable x:

y_{ij} = \alpha_j + \beta_j * x_i

where x is normally distributed:

x_i \sim N(\mu_x, \sigma_x)

The residual covariances are assumed to be zero between the y_j.

The most common constraint to ensure identification is to fix one of the \beta_j to 1, fix the location of the latent variable to 0, and estimate the variance of the latent variable. I have implemented this in Stan, and it works well:

```
data{
int<lower=0> N;
int<lower=0> J;
matrix[N, J] y;
}
parameters{
vector[J] alpha;
vector[J-1] beta;
vector<lower=0>[J] sigma;
vector[N] x;
real<lower=0> x_sigma;
}
transformed parameters{
// the fixed parameter
real beta_fix[J];
beta_fix[1] = 1;
for(j in 1:J-1){
beta_fix[j+1] = beta[j];
}
}
model{
matrix[N,J] mu;
x_sigma ~ normal(0, 1);
alpha ~ normal(0, 10);
beta ~ normal(0, 1);
sigma ~ normal(0, 1);
// latent variable
x ~ normal(0, x_sigma);
for(i in 1:N){
for(j in 1:J){
mu[i,j] = alpha[j] + beta_fix[j] * x[i];
y[i,j] ~ normal(mu[i,j], sigma[j]);
}
}
}
```

Running this model using some data from the `lavaan`

package using both lavaan and Stan:

```
require(lavaan)
data("HolzingerSwineford1939")
# lavaan model
HS.model <- 'visual =~ x1 + x2 + x3'
lavaan_fit1 <- cfa(HS.model, data=HolzingerSwineford1939, meanstructure=TRUE)
summary(lavaan_fit1)
# the Stan model
y <- as.matrix(HolzingerSwineford1939[,7:9])
N <- dim(y)[1]
J <- dim(y)[2]
stan_data <- list(
y = y, N = N, J = J
)
# fixed indicator/beta coefficient model
stan_fit_ff <- stan(file = "Stan-LV-fixed-indicator.stan",
data = stan_data)
print(stan_fit_ff, pars=c("alpha","beta","sigma","x_sigma"))
```

A different parameterization is fix both the latent variableâ€™s location to 0 and scale 1, and estimate all the indicator coefficients.

I have tried implementing this is Stan but the \beta_j coefficients very infrequently converge, with effective sample sizes usually of 2. Sometimes the chains converge and I get identical estimates to `lavaan`

. The model code and R code to run this model is below.

```
data{
int<lower=0> N;
int<lower=0> J;
matrix[N, J] y;
}
parameters{
vector[J] alpha;
vector[J] beta;
vector<lower=0>[J] sigma;
vector[N] x;
}
transformed parameters{
// no fixed parameters
}
model{
matrix[N,J] mu;
alpha ~ normal(0, 10);
beta ~ normal(0, 1);
sigma ~ normal(0, 1);
// latent variable
x ~ normal(0, 1);
for(i in 1:N){
for(j in 1:J){
mu[i,j] = alpha[j] + beta[j] * x[i];
y[i,j] ~ normal(mu[i,j], sigma[j]);
}
}
}
```

```
require(lavaan)
data("HolzingerSwineford1939")
# lavaan model
HS.model <- 'visual =~ x1 + x2 + x3'
lavaan_fit2 <- cfa(HS.model, data=HolzingerSwineford1939, meanstructure=TRUE, std.lv=TRUE)
summary(lavaan_fit2)
# the Stan model
y <- as.matrix(HolzingerSwineford1939[,7:9])
N <- dim(y)[1]
J <- dim(y)[2]
stan_data <- list(
y = y, N = N, J = J
)
# standardized latent variable
stan_fit_std <- stan(file = "Stan-LV-std.stan",
data = stan_data)
print(stan_fit_std, pars=c("alpha","beta","sigma"))
```

Does anyone know why the second model is not converging well, and what I could do about it?

Thanks for your help!