How to address warnings about low ESS and low BFMI?

Modeling
1

I am attempting to fit a multilevel mediation model with random intercepts using rstan. However, after a large number of iterations, I am still encountering warnings about low ESS and low BFMI. How should I address these warnings? Any response would be greatly appreciated.

My code goes here:
modelstanid11="data{
int<lower=1> N;
vector[N] c1;
vector[N] c2;
vector[N] aq;
vector[N] cr;
vector[N] adt;
int<lower=1> J;
int<lower=1, upper=J> Group[N];
int<lower=1> G;
int<lower=1, upper=G> Group2[N];
}
parameters{
real alphaaq;
real beta1aq;
real beta2aq;
real<lower=0> sigma_alphaaq;
real<lower=0> sigma_beta1aq;
real<lower=0> sigma_beta2aq;
real<lower=0> sigmaaq;

real alphacr;
real beta1cr;
real beta2cr;
real<lower=0> sigma_alphacr;
real<lower=0> sigma_beta1cr;
real<lower=0> sigma_beta2cr;
real<lower=0> sigmacr;

real alphaadt;
real beta1adt;
real beta2adt;
real beta3adt;
real beta4adt;
real<lower=0> sigma_alphaadt;
real<lower=0> sigma_beta1adt;
real<lower=0> sigma_beta2adt;
real<lower=0> sigma_beta3adt;
real<lower=0> sigma_beta4adt;
real<lower=0> sigma_adt;

cholesky_factor_corr[3] Omega_chol1;
cholesky_factor_corr[3] Omega_chol2;
vector<lower=0>[3] sigma_ranef1;
vector<lower=0>[3] sigma_ranef2;
matrix[J,3] gamma1; # random effects for id
matrix[G,3] gamma2; # random effects for topic
}
transformed parameters{
vector[J] gamma1alphaaq;

vector[J] gamma1alphacr;

vector[J] gamma1alpha;

for (j in 1:J){
gamma1alphaaq[j] = gamma1[j,1];

gamma1alphacr[j] = gamma1[j,2];

gamma1alpha[j]= gamma1[j,3];
}

vector[G] gamma2alphaaq;

vector[G] gamma2alphacr;

vector[G] gamma2alpha;

for (g in 1:G){
gamma2alphaaq[g] = gamma2[g,1];
gamma2alphacr[g] = gamma2[g,2];
gamma2alpha[g]= gamma2[g,3];}
}
model{
vector[N] mu_adt;
vector[N] mu_aq;
vector[N] mu_cr;
matrix[3,3] D1;
matrix[3,3] D2;
matrix[3, 3] DC1;
matrix[3, 3] DC2;

sigma_alphaaq ~ normal(0, 1);
sigma_beta1aq ~ normal(0, 1);
sigma_beta2aq ~ normal(0, 1);
alphaaq ~ normal(0, sigma_alphaaq);
beta1aq ~ normal(0, sigma_beta1aq);
beta2aq ~ normal(0, sigma_beta2aq);

sigmaaq ~ normal(0, 1);

sigma_alphacr ~ normal(0, 1);
sigma_beta1cr ~ normal(0, 1);
sigma_beta2cr ~ normal(0, 1);
alphacr ~ normal(0, sigma_alphacr);
beta1cr ~ normal(0, sigma_beta1cr);
beta2cr ~ normal(0, sigma_beta2cr);

sigmacr ~ normal(0, 1);

sigma_alphaadt ~ normal(0, 1);
sigma_beta1adt ~ normal(0, 1);
sigma_beta2adt ~ normal(0, 1);
sigma_beta3adt ~ normal(0, 1);
sigma_beta4adt ~ normal(0, 1);
alphaadt ~ normal(0, sigma_alphaadt);
beta1adt ~ normal(0, sigma_beta1adt);
beta2adt ~ normal(0, sigma_beta2adt);
beta3adt ~ normal(0, sigma_beta3adt);
beta4adt ~ normal(0, sigma_beta4adt);

sigma_adt ~ normal(0, 1);

sigma_ranef1 ~ normal(0, 1);
sigma_ranef2 ~ normal(0, 1);
Omega_chol1 ~ lkj_corr_cholesky(2.0);
Omega_chol2 ~ lkj_corr_cholesky(2.0);
D1 = diag_matrix(sigma_ranef1);
D2 = diag_matrix(sigma_ranef2);
DC1 = D1 * Omega_chol1;
DC2 = D2 * Omega_chol2;
for (j in 1:J) { # loop for Group random effects id
gamma1[j] ~ multi_normal_cholesky(rep_vector(0, 3), DC1);}
for (g in 1:G){ # loop for Group random effects topic
gamma2[g] ~ multi_normal_cholesky(rep_vector(0, 3), DC2);}

for (n in 1:N){
mu_aq[n] = alphaaq + gamma1alphaaq[Group[n]] + gamma2alphaaq[Group2[n]]+
(beta1aq) * c1[n] +
(beta2aq) * c2[n];

mu_cr[n] = alphacr + gamma1alphacr[Group[n]] + gamma2alphacr[Group2[n]]+
(beta1cr) * c1[n] +
(beta2cr) * c2[n];

mu_adt[n] = alphaadt + gamma1alpha[Group[n]] + gamma2alpha[Group2[n]]+
(beta1adt) * c1[n] +
(beta2adt) * c2[n] +
(beta3adt) * aq[n] +
(beta4adt) * cr[n];}

aq ~ normal(mu_aq, sigmaaq);
cr ~ normal(mu_cr, sigmacr);
adt ~ normal(mu_adt, sigma_adt);
}

generated quantities{
real naiveIndEffect1;

real naiveIndEffect2;

real naiveIndEffect3;

real naiveIndEffect4;

real totalEffect;

naiveIndEffect1=beta1aq*beta3adt;

naiveIndEffect2=beta1cr*beta4adt;

naiveIndEffect3=beta2aq*beta3adt;

naiveIndEffect4=beta2cr*beta4adt;

totalEffect=naiveIndEffect1+naiveIndEffect2+naiveIndEffect3+naiveIndEffect4+beta1adt
+beta2adt;
}
"

2

I’ve come across some cases, for instance, using reparameterization methods to address these issues, but I still don’t have a clear solution. :(

3

Unfortunately, there isn’t always a programatic way to fix issues such as these, and it often requires experimentation and tinkering. Good general purpose advice is presented here Runtime warnings and convergence problems

1 Like
5

Thank you very much for your response! I have read those general methods, but I still don’t know how to improve my model. The issue of low effective sample size (ess) has disappeared after a large number of iterations (60,000 times). Now, the most challenging problem is the low Bayesian Fraction of Missing Information (bfmi), and this warning still persists. I tried replacing the heavy-tailed Cauchy distribution with Student’s t-distribution and exponential distribution, but unfortunately, it didn’t work. Should I consider trying non-centered parameterization? I really need some more specific advice, as a beginner in rstan, and would greatly appreciate any help again!

6

Assuming your model sampling fit is stored in stanfit, the following
script in R gives you a summary of your variables and stores it in
modelstanid11.diagnosis. Which variables shows low n_eff
and/or high Rhat values?
Those are the ones you may need to focus on. Take a traceplot
and pairsplot.

options(max.print=1000000, width = 999)
sink("modelstanid11.diagnosis")
print(summary(stanfit))
sink()
7

Thx!! Perhaps these more specific diagnostic insights can assist me in optimizing my code : ) I’ll get right to running it!
I will update these results on the forum.

8

When I compare your model with the model in the Stan manual,
https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html
we see that the Stan-Manual model uses a mean in the group-level predictors.
BTW, you don’t need to run the model with 60000 iterations, 200/200 iterations should be enough
for diagnosis.

9

I’ll take an example of the initial fit result (using a Cauchy distribution) and I used 30,000 iterations with 15,000 warm-up iterations.
1.My Rstan Code
modelstan="data{
int<lower=1> N;
vector[N] c1;
vector[N] c2;
vector[N] aq;
vector[N] cr;
vector[N] adt;
int<lower=1> J;
int<lower=1, upper=J> Group[N];
int<lower=1> G;
int<lower=1, upper=G> Group2[N];
}
parameters{
real alphaaq;
real beta1aq;
real beta2aq;
real<lower=0> sigma_alphaaq;
real<lower=0> sigma_beta1aq;
real<lower=0> sigma_beta2aq;
real<lower=0> sigmaaq;

real alphacr;
real beta1cr;
real beta2cr;
real<lower=0> sigma_alphacr;
real<lower=0> sigma_beta1cr;
real<lower=0> sigma_beta2cr;
real<lower=0> sigmacr;

real alphaadt;
real beta1adt;
real beta2adt;
real beta3adt;
real beta4adt;
real<lower=0> sigma_alphaadt;
real<lower=0> sigma_beta1adt;
real<lower=0> sigma_beta2adt;
real<lower=0> sigma_beta3adt;
real<lower=0> sigma_beta4adt;
real<lower=0> sigma_adt;

cholesky_factor_corr[3] Omega_chol1;
cholesky_factor_corr[3] Omega_chol2;
vector<lower=0>[3] sigma_ranef1;
vector<lower=0>[3] sigma_ranef2;
matrix[J,3] gamma1; # random effects for id
matrix[G,3] gamma2; # random effects for topic
}
transformed parameters{
vector[J] gamma1alphaaq;

vector[J] gamma1alphacr;

vector[J] gamma1alpha;

for (j in 1:J){
gamma1alphaaq[j] = gamma1[j,1];

gamma1alphacr[j] = gamma1[j,2];

gamma1alpha[j]= gamma1[j,3];
}

vector[G] gamma2alphaaq;

vector[G] gamma2alphacr;

vector[G] gamma2alpha;

for (g in 1:G){
gamma2alphaaq[g] = gamma2[g,1];
gamma2alphacr[g] = gamma2[g,2];
gamma2alpha[g]= gamma2[g,3];}
}
model{
vector[N] mu_adt;
vector[N] mu_aq;
vector[N] mu_cr;
matrix[3,3] D1;
matrix[3,3] D2;
matrix[3, 3] DC1;
matrix[3, 3] DC2;

sigma_alphaaq ~ cauchy(0, 1);
sigma_beta1aq ~ cauchy(0, 1);
sigma_beta2aq ~ cauchy(0, 1);
alphaaq ~ normal(0, sigma_alphaaq);
beta1aq ~ normal(0, sigma_beta1aq);
beta2aq ~ normal(0, sigma_beta2aq);

sigmaaq ~ cauchy(0, 1);

sigma_alphacr ~ cauchy(0, 1);
sigma_beta1cr ~ cauchy(0, 1);
sigma_beta2cr ~ cauchy(0, 1);
alphacr ~ normal(0, sigma_alphacr);
beta1cr ~ normal(0, sigma_beta1cr);
beta2cr ~ normal(0, sigma_beta2cr);

sigmacr ~ cauchy(0, 1);

sigma_alphaadt ~ cauchy(0, 1);
sigma_beta1adt ~ cauchy(0, 1);
sigma_beta2adt ~ cauchy(0, 1);
sigma_beta3adt ~ cauchy(0, 1);
sigma_beta4adt ~ cauchy(0, 1);
alphaadt ~ normal(0, sigma_alphaadt);
beta1adt ~ normal(0, sigma_beta1adt);
beta2adt ~ normal(0, sigma_beta2adt);
beta3adt ~ normal(0, sigma_beta3adt);
beta4adt ~ normal(0, sigma_beta4adt);

sigma_adt ~ cauchy(0, 1);

sigma_ranef1 ~ cauchy(0, 1);
sigma_ranef2 ~ cauchy(0, 1);
Omega_chol1 ~ lkj_corr_cholesky(2.0);
Omega_chol2 ~ lkj_corr_cholesky(2.0);
D1 = diag_matrix(sigma_ranef1);
D2 = diag_matrix(sigma_ranef2);
DC1 = D1 * Omega_chol1;
DC2 = D2 * Omega_chol2;
for (j in 1:J) { # loop for Group random effects id
gamma1[j] ~ multi_normal_cholesky(rep_vector(0, 3), DC1);}
for (g in 1:G){ # loop for Group random effects topic
gamma2[g] ~ multi_normal_cholesky(rep_vector(0, 3), DC2);}

for (n in 1:N){
mu_aq[n] = alphaaq + gamma1alphaaq[Group[n]] + gamma2alphaaq[Group2[n]]+
(beta1aq) * c1[n] +
(beta2aq) * c2[n];

mu_cr[n] = alphacr + gamma1alphacr[Group[n]] + gamma2alphacr[Group2[n]]+
(beta1cr) * c1[n] +
(beta2cr) * c2[n];

mu_adt[n] = alphaadt + gamma1alpha[Group[n]] + gamma2alpha[Group2[n]]+
(beta1adt) * c1[n] +
(beta2adt) * c2[n] +
(beta3adt) * aq[n] +
(beta4adt) * cr[n];}

aq ~ normal(mu_aq, sigmaaq);
cr ~ normal(mu_cr, sigmacr);
adt ~ normal(mu_adt, sigma_adt);
}

generated quantities{
real naiveIndEffect1;

real naiveIndEffect2;

real naiveIndEffect3;

real naiveIndEffect4;

real totalEffect;

naiveIndEffect1=beta1aq*beta3adt;

naiveIndEffect2=beta1cr*beta4adt;

naiveIndEffect3=beta2aq*beta3adt;

naiveIndEffect4=beta2cr*beta4adt;

totalEffect=naiveIndEffect1+naiveIndEffect2+naiveIndEffect3+naiveIndEffect4+beta1adt
+beta2adt;
}
"
fit1= stan(model_code = modelstan, data=data_list10, iter = 400, warmup=200,
thin=1, chains = 1, verbose = F)
iter = 30000
wu = 15000
thin = 20
chains = 4
fit = stan(modelstan, fit=fit1,
data=data_list, iter=iter, warmup=wu,
thin=thin, cores=8,
control = list(adapt_delta = 0.9999, max_treedepth = 20))

2.The Code Used for Extracting Results
mainpars2=c(‘alphaaq’,
‘alphacr’,
‘alphaadt’,
‘beta1aq’,
‘beta2aq’,
‘beta1cr’,
‘beta2cr’,
‘alphaadt’,
‘beta1adt’,
‘beta2adt’,Code / Program output
‘beta3adt’,
‘beta4adt’,
‘sigmaaq’,
‘sigmacr’,
‘sigma_adt’,
‘sigma_ranef1’,
‘sigma_ranef2’,
‘naiveIndEffect1’,
‘naiveIndEffect2’,
‘naiveIndEffect3’,
‘naiveIndEffect4’,
‘totalEffect’)
pairs(fit, pars = c( “sigma_alphaaq”, “sigma_beta1aq”,“sigma_beta2aq”,“alphaaq”,“beta1aq”,“beta2aq”))
pairs(fit, pars = c( “sigma_alphacr”, “sigma_beta1cr”,“sigma_beta2cr”,“alphacr”,“beta1cr”,“beta2cr”))
pairs(fit, pars = c( “sigma_alphaadt”, “sigma_beta1adt”,“sigma_beta2adt”,“sigma_beta3adt”,“sigma_beta4adt”,“alphaadt”,“beta1adt”,“beta2adt”,“beta3adt”,“beta4adt”))
pairs(fit, pars = c( “sigma_ranef1”, “sigma_ranef2”))

print(fit,digits=3, probs = c(.025, .5, 0.975),pars=mainpars2)
3.Result
Inference for Stan model: anon_model.
4 chains, each with iter=30000; warmup=15000; thin=20;
post-warmup draws per chain=750, total post-warmup draws=3000.

              mean se_mean    sd   2.5%    50%  97.5% n_eff  Rhat

alphaaq -0.277 0.006 0.265 -0.803 -0.275 0.196 1800 1.003
alphacr -0.276 0.005 0.258 -0.764 -0.269 0.188 3055 0.999
alphaadt 0.036 0.001 0.051 -0.060 0.036 0.136 1567 1.002
beta1aq 0.398 0.001 0.041 0.319 0.398 0.477 1966 1.002
beta2aq 0.495 0.001 0.040 0.418 0.493 0.575 2826 1.000
beta1cr 0.378 0.002 0.041 0.298 0.378 0.456 356 1.016
beta2cr 0.473 0.001 0.040 0.393 0.473 0.548 1034 1.007
beta1adt -0.030 0.000 0.014 -0.057 -0.031 -0.002 2920 1.000
beta2adt -0.042 0.000 0.014 -0.070 -0.043 -0.015 2917 1.001
beta3adt 0.733 0.001 0.030 0.674 0.732 0.793 2824 1.001
beta4adt 0.237 0.001 0.030 0.176 0.238 0.296 2770 1.001
sigmaaq 0.746 0.001 0.015 0.720 0.746 0.776 646 1.007
sigmacr 0.743 0.000 0.015 0.715 0.742 0.772 2683 1.000
sigma_adt 0.243 0.000 0.005 0.232 0.243 0.253 889 1.005
sigma_ranef1[1] 0.444 0.000 0.027 0.392 0.444 0.500 3052 0.999
sigma_ranef1[2] 0.476 0.001 0.028 0.424 0.475 0.531 1038 1.004
sigma_ranef1[3] 0.094 0.000 0.010 0.073 0.095 0.114 2384 1.000
sigma_ranef2[1] 0.573 0.006 0.302 0.247 0.505 1.343 2549 1.000
sigma_ranef2[2] 0.556 0.006 0.297 0.233 0.476 1.313 2141 1.001
sigma_ranef2[3] 0.095 0.001 0.075 0.030 0.074 0.294 2526 1.000
naiveIndEffect1 0.292 0.001 0.032 0.231 0.292 0.356 2876 1.001
naiveIndEffect2 0.089 0.000 0.015 0.062 0.088 0.120 1720 1.005
naiveIndEffect3 0.363 0.001 0.033 0.301 0.361 0.429 2719 1.000
naiveIndEffect4 0.112 0.001 0.017 0.079 0.112 0.146 1104 1.005
totalEffect 0.784 0.001 0.048 0.690 0.785 0.875 2935 1.000

adt
adt1056×1042 167 KB

aq
aq977×1409 175 KB

cr
cr977×1409 163 KB

ranef
ranef1056×1042 174 KB

10

Those sigma_ look not good, why have you changed it to cauchy(0, 1)?
What happens, if you run the model by

fit = stan(modelstan, fit=fit1,
data=data_list, iter=500, warmup=500,, cores=8,
control = list(adapt_delta = 0.8, max_treedepth = 20))

Then I would get rid of “single” sigma’s, replace

sigma_alphaaq ~ cauchy(0, 1);
alphaaq ~ normal(0, sigma_alphaaq);

by something like
alphaaq ~ normal(0, 5);

Next, I would systematically build the model. Initially, without any groups, then adding one group at a time, followed by removing and introducing the other group. I would observe its behavior at each step. If the model performs well, I would proceed to include additional parameters.

11

The code I followed for the Bayesian Multilevel Mediation is from Bayesian Multilevel Mediation | Model Estimation by Example, and it seems there are some differences. Does this mean that I should follow the steps in the rstan manual instead of this one?

12

The Stan User Manual is a extraordinary good source. I almost build up my model from scratch, made tons of mistakes in the past and get experience what works and what not.
I am sure you may get the model to work. Just disable some of the feature and have a look how
it behaves.
This model itselfs could take some improvements, eg. the loop

for (n in 1:N){
mu_aq[n] = alphaaq + gamma1alphaaq[Group[n]] + gamma2alphaaq[Group2[n]]+
(beta1aq) * c1[n] +
(beta2aq) * c2[n];
}

can be improved by a vectorized version.

D2 = diag_matrix(sigma_ranef2);
DC1 = D1 * Omega_chol1;

The function diag_pre_multiply is usually a better choice.
The single sigma’s, unless you want shrinking. I don’t see, why one should
model it that way.
The following assignments are a unusual and not needed. One may use
gamma2 directly and matrix-multiplication.

for (g in 1:G){
gamma2alphaaq[g] = gamma2[g,1];
gamma2alphacr[g] = gamma2[g,2];
gamma2alpha[g]= gamma2[g,3];
}

But that doesn’t mean it’s an incorrect model, although it can be done better.

13

Great! These suggestions are extremely valuable, but I need time to understand and try them (which may not be an easy task for me). I will now proceed to make improvements step by step.

1 Like
15

data.csv (39.0 KB)
Here’s the data I used~

16

I read this case, in the section int<lower=1> L; // num group predictors , my case does not have group-level predictor variables (L=0), my model is similar to the following figure.

WechatIMG540
WechatIMG5401920×2560 213 KB

17

The following works without error, it uses the new Stan 3.x syntax. It should work with
rstan version 2.x using too, when you change back those two array statements.

  • Switch to the non-centered version
  • Introduced a mean parameter mu for the cov-distributions.
  • Use cmdstanr.

model_zrn_modern.R:

library(rstan)
library(cmdstanr)
fname <- "model_zrn_modern"
mod <- cmdstan_model(
  paste0(fname,".stan")
  , cpp_options = list(stan_threads = FALSE)
  , pedantic=T
  , force_recompile=F
  , compile_model_methods=TRUE 
  , stanc_options = list("O0")
)

csv_zrn <- read.table("data.csv",sep=",",header=T)
N <- max(csv_zrn$N)
c1 <- csv_zrn$c1
c2 <- csv_zrn$c2
J <- max(csv_zrn$Group)
G <- max(csv_zrn$Group2)
aq <- csv_zrn$aq
cr <- csv_zrn$cr
adt <- csv_zrn$adt

zrn_dat <- list(
    N      =  N
  , G      =  G
  , J      =  J
  , Group  =  csv_zrn$Group
  , Group2 =  csv_zrn$Group2
  , c1     =  c1
  , c2     =  c2
  , aq     =  aq
  , cr     =  cr
  , adt    =  adt
)

n_chains <- 2
###Sys.setenv(OPENBLAS_NUM_THREADS="1")
zrn_fit <- mod$sample(
    data = zrn_dat
  , seed = 12
  , chains = n_chains
  , parallel_chains = n_chains
  , refresh = 1
  , iter_warmup = 500
  , iter_sampling = 500
  , max_treedepth = 16
  , adapt_delta = 0.99
  # , step_size = 0.1
)
zrn_samp <- rstan::read_stan_csv(zrn_fit$output_files())
extr<-rstan::extract(zrn_samp)


options(max.print=1000000, width = 999)
sink(paste0(fname,".diagnosis"))
print(summary(zrn_fit))
sink()

model_zrn_modern.stan:

data{
  int<lower=1> N;
  vector[N] c1;
  vector[N] c2;
  vector[N] aq;
  vector[N] cr;
  vector[N] adt;
  int<lower=1> J;
  array[N] int<lower=1, upper=J> Group;
  int<lower=1> G;
  array[N] int<lower=1, upper=G> Group2;
}
parameters{
  real alphaaq;
  real beta1aq;
  real beta2aq;
  real<lower=0> sigma_alphaaq;
  real<lower=0> sigma_beta1aq;
  real<lower=0> sigma_beta2aq;
  real<lower=0> sigmaaq;
  real alphacr;
  real beta1cr;
  real beta2cr;
  real<lower=0> sigma_alphacr;
  real<lower=0> sigma_beta1cr;
  real<lower=0> sigma_beta2cr;
  real<lower=0> sigmacr;
  real alphaadt;
  real beta1adt;
  real beta2adt;
  real beta3adt;
  real beta4adt;
  real<lower=0> sigma_alphaadt;
  real<lower=0> sigma_beta1adt;
  real<lower=0> sigma_beta2adt;
  real<lower=0> sigma_beta3adt;
  real<lower=0> sigma_beta4adt;
  real<lower=0> sigma_adt;
  cholesky_factor_corr[3] Omega_chol1;
  cholesky_factor_corr[3] Omega_chol2;
  vector<lower=0>[3] sigma_ranef1;
  vector<lower=0>[3] sigma_ranef2;
  matrix[J,3] gamma1_raw;  // random effects for id
  matrix[G,3] gamma2_raw;  // random effects for topic
  vector[3] mu1;
  vector[3] mu2;
}
transformed parameters{
  vector[J] gamma1alphaaq;
  vector[J] gamma1alphacr;
  vector[J] gamma1alpha;
  vector[G] gamma2alphaaq;
  vector[G] gamma2alphacr;
  vector[G] gamma2alpha;
  vector[N] mu_adt;
  vector[N] mu_aq;
  vector[N] mu_cr;
  matrix[J,3] gamma1 = rep_matrix(mu1', J) + gamma1_raw * diag_pre_multiply(sigma_ranef1, Omega_chol1)';  // random effects for id
  matrix[G,3] gamma2 = rep_matrix(mu2', G) + gamma2_raw * diag_pre_multiply(sigma_ranef2, Omega_chol2)';

  // matrix[3,3] DC1 = diag_pre_multiply(sigma_ranef1, Omega_chol1);
  // matrix[3,3] DC2 = diag_pre_multiply(sigma_ranef2, Omega_chol2);
  for (j in 1:J){
    gamma1alphaaq[j] = gamma1[j,1];
    gamma1alphacr[j] = gamma1[j,2];
    gamma1alpha[j]= gamma1[j,3];
  }
  for (g in 1:G) {
    gamma2alphaaq[g] = gamma2[g,1];
    gamma2alphacr[g] = gamma2[g,2];
    gamma2alpha[g]= gamma2[g,3];
  }
  for (n in 1:N) {
    mu_aq[n] = alphaaq + gamma1alphaaq[Group[n]] + gamma2alphaaq[Group2[n]]+
    (beta1aq) * c1[n] +
    (beta2aq) * c2[n];
    mu_cr[n] = alphacr + gamma1alphacr[Group[n]] + gamma2alphacr[Group2[n]]+
    (beta1cr) * c1[n] +
    (beta2cr) * c2[n];
    mu_adt[n] = alphaadt + gamma1alpha[Group[n]] + gamma2alpha[Group2[n]]+
    (beta1adt) * c1[n] +
    (beta2adt) * c2[n] +
    (beta3adt) * aq[n] +
    (beta4adt) * cr[n];
  }
}
model{
  sigma_alphaaq ~ student_t(5, 0, 2);
  sigma_beta1aq ~ student_t(5, 0, 2);
  sigma_beta2aq ~ student_t(5, 0, 2);
  alphaaq ~ normal(0, sigma_alphaaq);
  beta1aq ~ normal(0, sigma_beta1aq);
  beta2aq ~ normal(0, sigma_beta2aq);
  sigmaaq ~ student_t(5, 0, 2);
  sigma_alphacr ~ student_t(5, 0, 2);
  sigma_beta1cr ~ student_t(5, 0, 2);
  sigma_beta2cr ~ student_t(5, 0, 2);
  alphacr ~ normal(0, sigma_alphacr);
  beta1cr ~ normal(0, sigma_beta1cr);
  beta2cr ~ normal(0, sigma_beta2cr);
  to_vector(gamma1_raw) ~ normal(0, 1);
  to_vector(gamma2_raw) ~ normal(0, 1);
  sigmacr ~ student_t(5, 0, 2);
  sigma_alphaadt ~ student_t(5, 0, 2);
  sigma_beta1adt ~ student_t(5, 0, 2);
  sigma_beta2adt ~ student_t(5, 0, 2);
  sigma_beta3adt ~ student_t(5, 0, 2);
  sigma_beta4adt ~ student_t(5, 0, 2);
  alphaadt ~ normal(0, sigma_alphaadt);
  beta1adt ~ normal(0, sigma_beta1adt);
  beta2adt ~ normal(0, sigma_beta2adt);
  beta3adt ~ normal(0, sigma_beta3adt);
  beta4adt ~ normal(0, sigma_beta4adt);
  sigma_adt ~ student_t(5, 0, 2);
  sigma_ranef1 ~ student_t(5, 0, 2);
  sigma_ranef2 ~ student_t(5, 0, 2);
  Omega_chol1 ~ lkj_corr_cholesky(2.0);
  Omega_chol2 ~ lkj_corr_cholesky(2.0);
  aq ~ normal(mu_aq, sigmaaq);
  cr ~ normal(mu_cr, sigmacr);
  adt ~ normal(mu_adt, sigma_adt);
  mu1 ~ normal(0, 2);
  mu2 ~ normal(0, 2);
  // for (j in 1:J) 
  //   gamma1[j] ~ multi_normal_cholesky(mu1, DC1);
  // for (g in 1:G)
  //   gamma2[g] ~ multi_normal_cholesky(mu2, DC2);
}

generated quantities{
  real naiveIndEffect1;
  real naiveIndEffect2;
  real naiveIndEffect3;
  real naiveIndEffect4;
  real totalEffect;
  naiveIndEffect1=beta1aq*beta3adt;
  naiveIndEffect2=beta1cr*beta4adt;
  naiveIndEffect3=beta2aq*beta3adt;
  naiveIndEffect4=beta2cr*beta4adt;
  totalEffect=naiveIndEffect1+naiveIndEffect2+naiveIndEffect3+naiveIndEffect4+beta1adt
  +beta2adt;
}
1 Like
19

Thank you so much for your help, you are very professional! I will check and run it now:)

20

Oh! I haven’t downloaded cmdstanr yet. I’m using the following code to install it (for macOS), but it seems to take some time.

21

if (!requireNamespace(“remotes”, quietly = TRUE))
install.packages(“remotes”)
remotes::install_github(“stan-dev/cmdstanr”)
install_cmdstan()

22

I should run with rstan, just change:

  array[N] int<lower=1, upper=J> Group;
  array[N] int<lower=1, upper=G> Group2;

to the original version.
I wanted to have a test with the newest Stan version 2.34.1.

24

1.I am using cmdstanr for analysis, but I encountered an error message:

mod ← cmdstan_model(

  • paste0(fname,“.stan”)
  • , cpp_options = list(stan_threads = FALSE)
  • , pedantic=T
  • , force_recompile=F
  • , compile_model_methods=TRUE
  • , stanc_options = list(“O0”)
  • )
    Compiling Stan program…
    Error in process_initialize(self, private, command, args, stdin, stdout, …:
    ! Native call to processx_exec failed
    Caused by error in chain_call(c_processx_exec, command, c(command, args), pty, pty_options, …:
    ! cannot start processx process ‘bin/stanc’ (system error 2, No such file or directory) @unix/processx.c:613 (processx_exec)
    Type .Last.error to see the more details.

This looks like an installation issue with cmdstan?

2.When I ran model_zrn_modern.stan using rstan, I changed

array[N] int<lower=1, upper=J> Group;
array[N] int<lower=1, upper=G> Group2;

to

int<lower=1, upper=J> Group[N];
int<lower=1, upper=G> Group2[N];

I got the warning:
"1: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat
2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
Runtime warnings and convergence problems "
THE CODE I USED:

test0126 = stan(model_code =model_zrn_modern, data=zrn_dat, iter = 30000, warmup=15000,
thin=1, chains = 1, verbose = F,control = list(max_treedepth = 20))

Does this mean my approach to running it is incorrect? Could you please provide me with the command you used to test this rstan model? It seems like our results are a bit inconsistent.

3. Is there an extra “#” symbol in this place?

# , step_size = 0.1

Thank you for your patience and all the help!