Hi all,

I am trying to fit a dynamic latent variable model with 2 latent variables that evolve over time. But I’m having some trouble getting the sampler to work (the trace plots look almost flat for several chains). Specifically, the sampler warns me that the E-BFMI is low (e.g. 0.00535) and that I may need to reparameterize my model. But I’ve already tried the reparametization based on the stochastic volatility model (i.e. https://mc-stan.org/docs/2_19/stan-users-guide/stochastic-volatility-models.html) and I’m simply not sure how to proceed from here.

I’ve also imposed sign constraints because I observed that there was sign indeterminacy for some of the parameters while trying different variants of the model.

The idea is that the second latent variable (i.e. F) should decay as a function of time while the first latent variable should describe the temporal evolution of the variables T_t, R_t, I_t . R_t is actually reaction time and its coded as -1 when there is no reaction. I_t is the indicator for reaction/no reaction.

Is the poor sampling behavior due to unidentifiability? Or is it because the two latent share some descendants in the graphical model (like how cross-loadings often lead to estimation issues in factor analysis)?

Would greatly appreciate any insights into this issue!

```
data {
int<lower=0> N;
int<lower=1> numBlocks;
int<lower=1> num_timepts;
int<lower=0,upper=1> I[numBlocks, N];
matrix<lower=-1>[numBlocks, N] R;
matrix<lower=0,upper=1>[numBlocks, N] V;
matrix<lower=0,upper=100>[numBlocks, N] T;
vector[N] G;
int<lower=0,upper=1> E[numBlocks, num_timepts, N];
int<lower=0,upper=1> S[numBlocks, num_timepts, N];
matrix[N,2] Grp;
matrix<lower=0>[numBlocks, N] Z;
}
transformed data {
int num_S0[numBlocks, N];
int num_S1[numBlocks, N];
int num_E_S0[numBlocks, N];
int num_E_S1[numBlocks, N];
matrix<lower=0>[numBlocks, N] I_saw_event;
for (blk_i in 1:numBlocks) {
for (per_i in 1:N) {
num_S0[blk_i, per_i] = 0;
num_S1[blk_i, per_i] = 0;
num_E_S0[blk_i, per_i] = 0;
num_E_S1[blk_i, per_i] = 0;
for (time_i in 1:num_timepts) {
num_S0[blk_i, per_i] += (1 - S[blk_i, time_i, per_i]);
if ( S[blk_i, time_i, per_i] == 0 ) {
num_E_S0[blk_i, per_i] += E[blk_i, time_i, per_i];
} else {
num_E_S1[blk_i, per_i] += E[blk_i, time_i, per_i];
}
}
num_S1[blk_i, per_i] = num_timepts - num_S0[blk_i, per_i];
}
}
I_saw_event = 1.0 - to_matrix(I);
}
parameters {
row_vector[N] B1_std;
row_vector[N] F1_std;
real b_B_Bt_slope;
real b_B_Tt_slope_raw; // Constrained to be positive and set upper bound
real b_B_Rt_slope_raw; // Constrained to be positive
real b_V_Bt_slope_raw; // Constrained to be positive
real b_V_It_slope;
real b_V_It_int;
real b_F_It_slope_raw; // Constrained to be negative
real b_F_Et_slope_raw; // Constrained positive
real F_int;
real b_Grp_Ft_slope_1;
real b_Grp_Ft_slope_2_raw; // Constrain to be negative
real b_Z_Bt_slope_raw; // Constrain to be positive
real<lower=0> sigma_Rt;
real<lower=0, upper=2> sigma_Tt; // Impose upper bound
real<lower=0> sigma_G;
real F_G_slope;
}
transformed parameters {
matrix[numBlocks+1, N] B;
matrix[numBlocks, N] F;
real<lower=5,upper=100> b_B_Tt_slope;
real<lower=0> b_B_Rt_slope;
real<lower=0> b_V_Bt_slope;
real<upper=0> b_F_It_slope;
vector<upper=0>[N] F_decay;
real<upper=0> b_Grp_Ft_slope_2;
row_vector[3] state_coefs;
real<lower=0> b_Z_Bt_slope;
real<lower=0> b_F_Et_slope;
b_B_Tt_slope = 95 * inv_logit(b_B_Tt_slope_raw) + 5;
b_B_Rt_slope = exp(b_B_Rt_slope_raw);
b_V_Bt_slope = exp(b_V_Bt_slope_raw);
b_Z_Bt_slope = exp(b_Z_Bt_slope_raw);
b_Grp_Ft_slope_2 = -exp(b_Grp_Ft_slope_2_raw);
b_F_It_slope = -exp(b_F_It_slope_raw);
b_F_Et_slope = exp(b_F_Et_slope_raw);
for (per_i in 1:N) {
F_decay[per_i] = -exp(b_Grp_Ft_slope_1 * (Grp'[1])'[per_i] + b_Grp_Ft_slope_2 * (Grp'[2])'[per_i] + F_int);
}
// Initialize
B[1] = B1_std;
F[1] = exp(F1_std);
// State equation
state_coefs = [b_V_Bt_slope, b_B_Bt_slope, b_Z_Bt_slope];
for (t in 2:(numBlocks+1)) {
B[t] = state_coefs * [V[t-1] .* I_saw_event[t-1], B[t-1], Z[t-1]];
if (t <= numBlocks) {
for (per_i in 1:N) {
F[t, per_i] = pow(F[1, per_i], (t-1) * F_decay[per_i]);
}
}
}
}
model {
// Priors for coefficients
F_int ~ normal(0,5);
b_Grp_Ft_slope_1 ~ normal(0,5);
b_Grp_Ft_slope_2 ~ normal(0,5);
F_G_slope ~ normal(0,5);
b_B_Bt_slope ~ normal(0,5);
b_V_It_slope ~ normal(0,5);
b_V_It_int ~ normal(0,5);
b_F_Et_slope_raw ~ normal(0,1);
b_F_It_slope_raw ~ normal(0,1);
sigma_Rt ~ normal(0,1);
sigma_Tt ~ normal(0,1);
b_B_Tt_slope_raw ~ normal(0,1);
b_B_Rt_slope_raw ~ normal(0,1);
b_V_Bt_slope_raw ~ normal(0,1);
b_Z_Bt_slope ~ normal(0,5);
sigma_G ~ normal(0,1);
F_G_slope ~ normal(0,5);
// Priors for latent variables
B1_std ~ std_normal();
F1_std ~ std_normal();
// Measurement model
for ( k in 1:numBlocks ) {
// Measurement for I
I[k] ~ bernoulli_logit( b_V_It_int + b_V_It_slope * V[k] + b_F_It_slope * F[k] );
// Measurement for R
for ( rt_n in 1:N ) {
if ( I[k, rt_n] == 0 ) {
R[k, rt_n] ~ lognormal( b_B_Rt_slope * B[k, rt_n] , sigma_Rt );
} else {
R[k, rt_n] ~ normal( -1, 0.001 );
}
}
// Measurement for T
T[k] ~ normal( b_B_Tt_slope * B[k+1], sigma_Tt );
// Measurement for E
for (e_n in 1:N) {
if ( num_S0[k, e_n] != 0 ) {
num_E_S0[k, e_n] ~ binomial_logit( num_S0[k, e_n], b_F_Et_slope * F[k, e_n] );
}
if ( num_S1[k, e_n] != 0 ) {
num_E_S1[k, e_n] ~ binomial_logit( num_S1[k, e_n], b_F_Et_slope * F[k, e_n] );
}
}
}
// Measurement for G
G ~ normal( F_G_slope * F' * [1./8., 1./8., 1./8., 1./8., 1./8., 1./8., 1./8., 1./8.]', sigma_G );
}
```