Divregeance appearing when incrasing chain number

Hello!
I am fitting a multivariate model with shared parameters. I am new to stan and had help from this forum to edit the stan code produced by brms (https://discourse.mc-stan.org/t/fit-multivariate-non-linear-model-with-shared-parameters/36243). I was able to fit the model relatively easily with two chains, thin = 1, adapt_delta = 0.9, and 8000 iterations. The fit was also very good. However, now that my model is performing well I want to run it with four chains to have a more final version. When I increase the number of chains to 3 or 4 I get wasnings about the low ESS, very high Rhat and very low n_eff. I tried stronger priors, more iteration and higher thin. It did not change anything. More iterations did not increase n_eff of the problemeatic parameters at all.

My stan code is:

functions {
}
data {
  int<lower=1> N;  // total number of observations
  int<lower=1> N_Nmr;  // number of observations
  vector[N_Nmr] Y_Nmr;  // response variable
  vector<lower=0>[N_Nmr] se_Nmr;  // known sampling error
  int<lower=1> K_Nmr_logI0;  // number of population-level effects
  matrix[N_Nmr, K_Nmr_logI0] X_Nmr_logI0;  // population-level design matrix
  int<lower=1> K_Nmr_b;  // number of population-level effects
  matrix[N_Nmr, K_Nmr_b] X_Nmr_b;  // population-level design matrix
  int<lower=1> K_Nmr_k;  // number of population-level effects
  matrix[N_Nmr, K_Nmr_k] X_Nmr_k;  // population-level design matrix
  int<lower=1> K_Nmr_int;  // number of population-level effects
  matrix[N_Nmr, K_Nmr_int] X_Nmr_int;  // population-level design matrix
  int<lower=1> K_Nmr_a;  // number of population-level effects
  matrix[N_Nmr, K_Nmr_a] X_Nmr_a;  // population-level design matrix
  // covariates for non-linear functions
  vector[N_Nmr] C_Nmr_1;
  vector[N_Nmr] C_Nmr_2;
  vector[N_Nmr] C_Nmr_3;
  vector[N_Nmr] C_Nmr_4;
  vector[N_Nmr] C_Nmr_5;
  vector[N_Nmr] C_Nmr_6;
  // covariates for non-linear functions
  vector[N_Nmr] C_sigma_Nmr_1;
  int<lower=1> N_BI;  // number of observations
  vector[N_BI] Y_BI;  // response variable
  vector<lower=0>[N_BI] se_BI;  // known sampling error
  int<lower=1> K_BI_logI0;  // number of population-level effects
  matrix[N_BI, K_BI_logI0] X_BI_logI0;  // population-level design matrix
  int<lower=1> K_BI_b;  // number of population-level effects
  matrix[N_BI, K_BI_b] X_BI_b;  // population-level design matrix
  int<lower=1> K_BI_k;  // number of population-level effects
  matrix[N_BI, K_BI_k] X_BI_k;  // population-level design matrix
  int<lower=1> K_BI_int;  // number of population-level effects
  matrix[N_BI, K_BI_int] X_BI_int;  // population-level design matrix
  int<lower=1> K_BI_a;  // number of population-level effects
  matrix[N_BI, K_BI_a] X_BI_a;  // population-level design matrix
  // covariates for non-linear functions
  vector[N_BI] C_BI_1;
  vector[N_BI] C_BI_2;
  vector[N_BI] C_BI_3;
  vector[N_BI] C_BI_4;
  vector[N_BI] C_BI_5;
  vector[N_BI] C_BI_6;
  // covariates for non-linear functions
  vector[N_BI] C_sigma_BI_1;
  int<lower=1> nresp;  // number of responses
  int nrescor;  // number of residual correlations
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  array[N_Nmr] int<lower=1> J_1_Nmr;  // grouping indicator per observation
  // group-level predictor values
  vector[N_Nmr] Z_1_Nmr_k_1;
  // data for group-level effects of ID 2
  int<lower=1> N_2;  // number of grouping levels
  int<lower=1> M_2;  // number of coefficients per level
  array[N_BI] int<lower=1> J_2_BI;  // grouping indicator per observation
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  vector<lower=0>[N_Nmr] se2_Nmr = square(se_Nmr);
  vector<lower=0>[N_BI] se2_BI = square(se_BI);
  array[N] vector[nresp] Y;  // response array
  for (n in 1:N) {
    Y[n] = transpose([Y_Nmr[n], Y_BI[n]]);
  }
}
parameters {
  vector<lower=0>[K_Nmr_logI0] b_Nmr_logI0;  // regression coefficients
  vector<lower=0,upper=1>[K_Nmr_b] b_Nmr_b;  // regression coefficients
  vector<lower=0>[K_Nmr_k] b_Nmr_k;  // regression coefficients
  vector<lower=0>[K_Nmr_int] b_Nmr_int;  // regression coefficients
  vector<lower=0>[K_Nmr_a] b_Nmr_a;  // regression coefficients
  vector<lower=0>[K_BI_int] b_BI_int;  // regression coefficients
  vector<lower=0>[K_BI_a] b_BI_a;  // regression coefficients
  cholesky_factor_corr[nresp] Lrescor;  // parameters for multivariate linear models
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  array[M_1] vector[N_1] z_1;  // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_Nmr_k_1;  // actual group-level effects
  real lprior = 0;  // prior contributions to the log posterior
  r_1_Nmr_k_1 = (sd_1[1] * (z_1[1]));
  lprior += normal_lpdf(b_Nmr_logI0 | 10, 2)
    - 1 * normal_lccdf(0 | 10, 2);
  lprior += beta_lpdf(b_Nmr_b | 0.7 * 4,(1 - 0.7) * 4);
  lprior += normal_lpdf(b_Nmr_k | 0,1)
    - 1 * normal_lccdf(0 | 0, 1);
  lprior += normal_lpdf(b_Nmr_int | 0, 1)
    - 1 * normal_lccdf(0 | 0, 1);
  lprior += normal_lpdf(b_Nmr_a | 0, 0.5)
    - 1 * normal_lccdf(0 | 0, 0.5);
  lprior += normal_lpdf(b_BI_int | 0, 1)
    - 1 * normal_lccdf(0 | 0, 1);
  lprior += normal_lpdf(b_BI_a | 0, 0.5)
    - 1 * normal_lccdf(0 | 0, 0.5);
  lprior += lkj_corr_cholesky_lpdf(Lrescor | 1);
  lprior += exponential_lpdf(sd_1 | 3);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N_Nmr] nlp_Nmr_logI0 = rep_vector(0.0, N_Nmr);
    // initialize linear predictor term
    vector[N_Nmr] nlp_Nmr_b = rep_vector(0.0, N_Nmr);
    // initialize linear predictor term
    vector[N_Nmr] nlp_Nmr_k = rep_vector(0.0, N_Nmr);
    // initialize linear predictor term
    vector[N_Nmr] nlp_Nmr_int = rep_vector(0.0, N_Nmr);
    // initialize linear predictor term
    vector[N_Nmr] nlp_Nmr_a = rep_vector(0.0, N_Nmr);
    // initialize non-linear predictor term
    vector[N_Nmr] mu_Nmr;
    // initialize non-linear predictor term
    vector[N_Nmr] sigma_Nmr;
    // initialize linear predictor term
    vector[N_BI] nlp_BI_int = rep_vector(0.0, N_BI);
    // initialize linear predictor term
    vector[N_BI] nlp_BI_a = rep_vector(0.0, N_BI);
    // initialize non-linear predictor term
    vector[N_BI] mu_BI;
    // initialize non-linear predictor term
    vector[N_BI] sigma_BI;
    // multivariate predictor array
    array[N] vector[nresp] Mu;
    array[N] vector[nresp] sigma;
    // cholesky factor of residual covariance matrix
    array[N] matrix[nresp, nresp] LSigma;
    nlp_Nmr_logI0 += X_Nmr_logI0 * b_Nmr_logI0;
    nlp_Nmr_b += X_Nmr_b * b_Nmr_b;
    nlp_Nmr_k += X_Nmr_k * b_Nmr_k;
    nlp_Nmr_int += X_Nmr_int * b_Nmr_int;
    nlp_Nmr_a += X_Nmr_a * b_Nmr_a;
    nlp_BI_int += X_BI_int * b_BI_int;
    nlp_BI_a += X_BI_a * b_BI_a;
    for (n in 1:N_Nmr) {
      // add more terms to the linear predictor
      nlp_Nmr_k[n] += r_1_Nmr_k_1[J_1_Nmr[n]] * Z_1_Nmr_k_1[n];
    }
    for (n in 1:N_Nmr) {
      // compute non-linear predictor values
      mu_Nmr[n] = ((C_Nmr_1[n] * (C_Nmr_2[n] + (inv_logit(nlp_Nmr_k[n]) * C_Nmr_3[n]))) / (exp(nlp_Nmr_logI0[n]) * C_Nmr_4[n] * (C_Nmr_5[n] ^ nlp_Nmr_b[n] + (C_Nmr_6[n] / 2) * nlp_Nmr_b[n] * (nlp_Nmr_b[n] - 1) * C_Nmr_5[n] ^ (nlp_Nmr_b[n] - 2))));
    }
    for (n in 1:N_Nmr) {
      // compute non-linear predictor values
      sigma_Nmr[n] = exp(nlp_Nmr_int[n] + C_sigma_Nmr_1[n] ^ nlp_Nmr_a[n]);
    }
    for (n in 1:N_BI) {
      // compute non-linear predictor values
      mu_BI[n] = (((C_BI_1[n] + inv_logit(nlp_Nmr_k[n]) * C_BI_2[n]) * C_BI_3[n] * C_BI_4[n]) / (exp(nlp_Nmr_logI0[n]) * C_BI_5[n] * (C_BI_4[n] ^ nlp_Nmr_b[n] + (C_BI_6[n] / 2) * nlp_Nmr_b[n] * (nlp_Nmr_b[n] - 1) * C_BI_4[n] ^ (nlp_Nmr_b[n] - 2))));
    }
    for (n in 1:N_BI) {
      // compute non-linear predictor values
      sigma_BI[n] = exp(nlp_BI_int[n] + C_sigma_BI_1[n] ^ nlp_BI_a[n]);
    }
    // combine univariate parameters
    for (n in 1:N) {
      Mu[n] = transpose([mu_Nmr[n], mu_BI[n]]);
      sigma[n] = transpose([sqrt(square(sigma_Nmr[n]) + se2_Nmr[n]), sqrt(square(sigma_BI[n]) + se2_BI[n])]);
      LSigma[n] = diag_pre_multiply(sigma[n], Lrescor);
    }
    for (n in 1:N) {
      target += multi_normal_cholesky_lpdf(Y[n] | Mu[n], LSigma[n]);
    }
  }
  // priors including constants
  target += lprior;
  target += std_normal_lpdf(z_1[1]);
}
generated quantities {
  // residual correlations
  corr_matrix[nresp] Rescor = multiply_lower_tri_self_transpose(Lrescor);
  vector<lower=-1,upper=1>[nrescor] rescor;
  // extract upper diagonal of correlation matrix
  for (k in 1:nresp) {
    for (j in 1:(k - 1)) {
      rescor[choose(k - 1, 2) + j] = Rescor[j, k];
    }
  }
}

Nmr and BI are highly correlated. It is a problem?

The stan data:

structure(list(N_Nmr = 35L, Y_Nmr = structure(c(145L, 41L, 72L, 
3173L, 2966L, 1262L, 1149L, 504L, 453L, 777L, 704L, 1923L, 1751L, 
2052L, 1921L, 195L, 167L, 245L, 208L, 805L, 327L, 456L, 78L, 
1747L, 936L, 934L, 326L, 1319L, 383L, 1350L, 82L, 1533L, 372L, 
896L, 193L), dim = 35L), se_Nmr = structure(c(36.4795918367347, 
33.9285714285714, 18.1122448979592, 336.989795918367, 336.989795918367, 
167.602040816327, 167.602040816327, 87.7551020408163, 87.7551020408163, 
146.428571428571, 146.428571428571, 233.163265306122, 233.163265306122, 
356.632653061225, 356.632653061225, 28.8265306122449, 28.8265306122449, 
38.7755102040816, 38.7755102040816, 70.9183673469388, 70.9183673469388, 
43.3673469387755, 43.3673469387755, 266.326530612245, 266.326530612245, 
116.326530612245, 116.326530612245, 109.948979591837, 109.948979591837, 
101.530612244898, 79.0816326530612, 163.520408163265, 163.520408163265, 
139.030612244898, 131.377551020408), dim = 35L), N_BI = 35L, 
    Y_BI = structure(c(58699.2067307692, 16338.5, 37062.72, 586326.283422458, 
    547110.145417652, 310087.009803921, 281911.414847162, 34735.4823529412, 
    31263.0941704036, 70328.0794520548, 63944.1487839771, 215919.163157895, 
    196235.196078432, 229360.696875, 214911.875, 26769.8888888888, 
    22947.037037037, 32365.8243243242, 27198.3756345178, 114220.875, 
    52637.9166666666, 42936.2584615384, 9754.75609756097, 75333.7133028107, 
    41971.8079673136, 45241.4505050505, 18686.2824207493, 126776.251643957, 
    44161.8114602588, 95337.7232142857, 6722.56559766764, 77728.6077844312, 
    21346.2576687117, 55276.7228915663, 9130.71748878923), dim = 35L), 
    se_BI = structure(c(14767.7455357143, 13520.5357142857, 9323.45918367347, 
    62271.0288660917, 62161.3406099672, 41181.6289015605, 41121.7828179307, 
    6048.04721888756, 6056.28260272719, 13253.5909980431, 13300.0715307582, 
    26180.1441102757, 26130.6905114988, 39862.3361766582, 39898.2780612245, 
    3957.34882842024, 3960.97883597882, 5122.45449531162, 5070.34082668601, 
    10062.5564868805, 11415.8871882086, 4083.40266875981, 5423.56271777003, 
    11484.4685186463, 11942.5277771987, 5634.66914038342, 6667.82332529554, 
    10567.7934076808, 12677.6660756724, 7170.14622813411, 6483.31052537633, 
    8291.07219235, 9383.19560953216, 8577.18373493977, 6215.39535096549
    ), dim = 35L), nresp = 2L, nrescor = 1, N = 35L, C_Nmr_1 = structure(c(592.2, 
    551.335652173913, 1720.40869551196, 5135.1, 3710.06869230388, 
    1904.81995681897, 1681.86864942529, 2445.9, 460.748600012, 
    2357.455533366, 1091.59457333333, 1240.4, 2162.041408085, 
    1302.2663333, 410.54, 3050.5, 2351.75561511765, 6834.90718954248, 
    5222.50784313725, 797.5, 847.063260001027, 696.263082154795, 
    1297.35043378995, 7805.1, 7483.01666657353, 12660.8274504265, 
    6406.01764705882, 917.4, 1138.47565652879, 1740.51025227273, 
    2371.19474747475, 1530.6, 4231.82828278788, 1926.86065651591, 
    619.783535353535), dim = 35L), C_Nmr_2 = structure(c(10662.3716311457, 
    1924.04403258699, 8519.96312725198, 44606.1973266835, 7148.80605247308, 
    40475.0499899619, 14328.9779815932, 1179.41546274407, 311.338271604938, 
    1011.65268657152, 343.723340046114, 30416.3327204273, 6934.1645973791, 
    28020.7274319236, 8795.01318623588, 15102.014475937, 2992.42688425335, 
    13679.0305628833, 4596.94257937868, 1431.80202918211, 321.916668837846, 
    1258.3342989412, 365.206573098678, 1384.29643593022, 221.813518739632, 
    1101.88876890077, 458.025245406228, 27277.3732597846, 6227.95201821207, 
    25551.6259521933, 6779.72156553717, 3858.32694159881, 849.379914193153, 
    3493.00224732218, 1107.93016373347), dim = 35L), C_Nmr_3 = structure(c(107809.545032839, 
    139852.603250838, 107809.545032839, 319346.983077122, 319346.983077122, 
    319346.983077122, 319346.983077122, 56377.3469017162, 142137.4817856, 
    142137.4817856, 142137.4817856, 438588.370047914, 438588.370047914, 
    438588.370047914, 1851206.799571, 19513.2105459205, 19513.2105459205, 
    19513.2105459205, 19513.2105459205, 320208.980098382, 320208.980098382, 
    320208.980098382, 320208.980098382, 37654.7145510446, 37654.7145510446, 
    37654.7145510446, 37654.7145510446, 355618, 330896.902438548, 
    330896.902438548, 330896.902438548, 156332.500004466, 60228.1080636982, 
    60228.1080636982, 60228.1080636982), dim = 35L), C_Nmr_4 = structure(c(0.968898826010513, 
    0.919258766903425, 0.952629817554708, 0.964098761243345, 
    0.915984845013819, 0.957674570493223, 0.900846934099921, 
    0.861315579722472, 0.891627183842253, 0.750831703960004, 
    0.900846934099921, 0.826518486768192, 0.640759038748287, 
    0.763040264421148, 0.774392686186511, 0.731328502676302, 
    0.68537658471524, 0.637465208746466, 0.68537658471524, 0.94955451815625, 
    0.80289108787463, 0.909417409482917, 0.794101192295099, 0.920525161523524, 
    0.677939358232424, 0.856408047245123, 0.76414932341953, 0.800447272830747, 
    0.626426585425339, 0.742450488957026, 0.620807014871961, 
    0.48341486108144, 0.624795727445052, 0.502013926753242, 0.620807014871961
    ), dim = 35L), C_Nmr_5 = structure(c(404.822115384615, 398.5, 
    514.76, 184.786096256684, 184.460601961447, 245.710784313725, 
    245.353711790393, 68.9196078431373, 69.0134529147982, 90.5123287671233, 
    90.829756795422, 112.282456140351, 112.070357554787, 111.77421875, 
    111.875, 137.281481481481, 137.407407407407, 132.105405405405, 
    130.761421319797, 141.889285714286, 160.972222222222, 94.1584615384615, 
    125.060975609756, 43.1217591887869, 44.8416751787538, 48.4383838383838, 
    57.3198847262248, 96.1154296011807, 115.304990757856, 70.6205357142857, 
    81.9825072886297, 50.7035928143713, 57.3824130879346, 61.6927710843374, 
    47.3094170403587), dim = 35L), C_Nmr_6 = structure(c(17205.5561081031, 
    7855.52631578947, 17981.9402222222, 7167.07851532402, 7083.85137352115, 
    4429.97998155698, 4310.38614224204, 1366.25000784314, 1255.20431299441, 
    563.231373668189, 586.130001516698, 5353.76039628483, 5295.71114620791, 
    2950.8776765502, 2976.15945805107, 5062.27618233618, 4841.68391994479, 
    6727.54052552553, 6431.56013674505, 19279.9067512909, 24098.3565351895, 
    14572.9452788462, 31025.1496333982, 505.623552253161, 431.953435227717, 
    464.226674912389, 321.986973397078, 26678.425812748, 53513.8419935647, 
    2334.9340789897, 3756.29209077115, 404.951071351273, 411.166985148681, 
    661.822873934763, 404.439865874843), dim = 35L), C_sigma_Nmr_1 = structure(c(592.2, 
    551.335652173913, 1720.40869551196, 5135.1, 3710.06869230388, 
    1904.81995681897, 1681.86864942529, 2445.9, 460.748600012, 
    2357.455533366, 1091.59457333333, 1240.4, 2162.041408085, 
    1302.2663333, 410.54, 3050.5, 2351.75561511765, 6834.90718954248, 
    5222.50784313725, 797.5, 847.063260001027, 696.263082154795, 
    1297.35043378995, 7805.1, 7483.01666657353, 12660.8274504265, 
    6406.01764705882, 917.4, 1138.47565652879, 1740.51025227273, 
    2371.19474747475, 1530.6, 4231.82828278788, 1926.86065651591, 
    619.783535353535), dim = 35L), K_Nmr_logI0 = 1L, X_Nmr_logI0 = structure(c(1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 
    1L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
    "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
    "28", "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_Nmr_b = 1L, X_Nmr_b = structure(c(1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_Nmr_k = 1L, X_Nmr_k = structure(c(1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    Z_1_Nmr_k_1 = structure(c(`1` = 1, `2` = 1, `3` = 1, `4` = 1, 
    `5` = 1, `6` = 1, `7` = 1, `8` = 1, `9` = 1, `10` = 1, `11` = 1, 
    `12` = 1, `13` = 1, `14` = 1, `15` = 1, `16` = 1, `17` = 1, 
    `18` = 1, `19` = 1, `20` = 1, `21` = 1, `22` = 1, `23` = 1, 
    `24` = 1, `25` = 1, `26` = 1, `27` = 1, `28` = 1, `29` = 1, 
    `30` = 1, `31` = 1, `32` = 1, `33` = 1, `34` = 1, `35` = 1
    ), dim = 35L, dimnames = list(c("1", "2", "3", "4", "5", 
    "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", 
    "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
    "27", "28", "29", "30", "31", "32", "33", "34", "35"))), 
    K_Nmr_int = 1L, X_Nmr_int = structure(c(1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_Nmr_a = 1L, X_Nmr_a = structure(c(1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    J_1_Nmr = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 
    7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L), dim = 35L), 
    C_BI_1 = structure(c(10662.3716311457, 1924.04403258699, 
    8519.96312725198, 44606.1973266835, 7148.80605247308, 40475.0499899619, 
    14328.9779815932, 1179.41546274407, 311.338271604938, 1011.65268657152, 
    343.723340046114, 30416.3327204273, 6934.1645973791, 28020.7274319236, 
    8795.01318623588, 15102.014475937, 2992.42688425335, 13679.0305628833, 
    4596.94257937868, 1431.80202918211, 321.916668837846, 1258.3342989412, 
    365.206573098678, 1384.29643593022, 221.813518739632, 1101.88876890077, 
    458.025245406228, 27277.3732597846, 6227.95201821207, 25551.6259521933, 
    6779.72156553717, 3858.32694159881, 849.379914193153, 3493.00224732218, 
    1107.93016373347), dim = 35L), C_BI_2 = structure(c(107809.545032839, 
    139852.603250838, 107809.545032839, 319346.983077122, 319346.983077122, 
    319346.983077122, 319346.983077122, 56377.3469017162, 142137.4817856, 
    142137.4817856, 142137.4817856, 438588.370047914, 438588.370047914, 
    438588.370047914, 1851206.799571, 19513.2105459205, 19513.2105459205, 
    19513.2105459205, 19513.2105459205, 320208.980098382, 320208.980098382, 
    320208.980098382, 320208.980098382, 37654.7145510446, 37654.7145510446, 
    37654.7145510446, 37654.7145510446, 355618, 330896.902438548, 
    330896.902438548, 330896.902438548, 156332.500004466, 60228.1080636982, 
    60228.1080636982, 60228.1080636982), dim = 35L), C_BI_3 = structure(c(592.2, 
    551.335652173913, 1720.40869551196, 5135.1, 3710.06869230388, 
    1904.81995681897, 1681.86864942529, 2445.9, 460.748600012, 
    2357.455533366, 1091.59457333333, 1240.4, 2162.041408085, 
    1302.2663333, 410.54, 3050.5, 2351.75561511765, 6834.90718954248, 
    5222.50784313725, 797.5, 847.063260001027, 696.263082154795, 
    1297.35043378995, 7805.1, 7483.01666657353, 12660.8274504265, 
    6406.01764705882, 917.4, 1138.47565652879, 1740.51025227273, 
    2371.19474747475, 1530.6, 4231.82828278788, 1926.86065651591, 
    619.783535353535), dim = 35L), C_BI_4 = structure(c(404.822115384615, 
    398.5, 514.76, 184.786096256684, 184.460601961447, 245.710784313725, 
    245.353711790393, 68.9196078431373, 69.0134529147982, 90.5123287671233, 
    90.829756795422, 112.282456140351, 112.070357554787, 111.77421875, 
    111.875, 137.281481481481, 137.407407407407, 132.105405405405, 
    130.761421319797, 141.889285714286, 160.972222222222, 94.1584615384615, 
    125.060975609756, 43.1217591887869, 44.8416751787538, 48.4383838383838, 
    57.3198847262248, 96.1154296011807, 115.304990757856, 70.6205357142857, 
    81.9825072886297, 50.7035928143713, 57.3824130879346, 61.6927710843374, 
    47.3094170403587), dim = 35L), C_BI_5 = structure(c(0.968898826010513, 
    0.919258766903425, 0.952629817554708, 0.964098761243345, 
    0.915984845013819, 0.957674570493223, 0.900846934099921, 
    0.861315579722472, 0.891627183842253, 0.750831703960004, 
    0.900846934099921, 0.826518486768192, 0.640759038748287, 
    0.763040264421148, 0.774392686186511, 0.731328502676302, 
    0.68537658471524, 0.637465208746466, 0.68537658471524, 0.94955451815625, 
    0.80289108787463, 0.909417409482917, 0.794101192295099, 0.920525161523524, 
    0.677939358232424, 0.856408047245123, 0.76414932341953, 0.800447272830747, 
    0.626426585425339, 0.742450488957026, 0.620807014871961, 
    0.48341486108144, 0.624795727445052, 0.502013926753242, 0.620807014871961
    ), dim = 35L), C_BI_6 = structure(c(17205.5561081031, 7855.52631578947, 
    17981.9402222222, 7167.07851532402, 7083.85137352115, 4429.97998155698, 
    4310.38614224204, 1366.25000784314, 1255.20431299441, 563.231373668189, 
    586.130001516698, 5353.76039628483, 5295.71114620791, 2950.8776765502, 
    2976.15945805107, 5062.27618233618, 4841.68391994479, 6727.54052552553, 
    6431.56013674505, 19279.9067512909, 24098.3565351895, 14572.9452788462, 
    31025.1496333982, 505.623552253161, 431.953435227717, 464.226674912389, 
    321.986973397078, 26678.425812748, 53513.8419935647, 2334.9340789897, 
    3756.29209077115, 404.951071351273, 411.166985148681, 661.822873934763, 
    404.439865874843), dim = 35L), C_sigma_BI_1 = structure(c(592.2, 
    551.335652173913, 1720.40869551196, 5135.1, 3710.06869230388, 
    1904.81995681897, 1681.86864942529, 2445.9, 460.748600012, 
    2357.455533366, 1091.59457333333, 1240.4, 2162.041408085, 
    1302.2663333, 410.54, 3050.5, 2351.75561511765, 6834.90718954248, 
    5222.50784313725, 797.5, 847.063260001027, 696.263082154795, 
    1297.35043378995, 7805.1, 7483.01666657353, 12660.8274504265, 
    6406.01764705882, 917.4, 1138.47565652879, 1740.51025227273, 
    2371.19474747475, 1530.6, 4231.82828278788, 1926.86065651591, 
    619.783535353535), dim = 35L), K_BI_logI0 = 1L, X_BI_logI0 = structure(c(1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 
    1L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
    "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
    "28", "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_BI_b = 1L, X_BI_b = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_BI_k = 1L, X_BI_k = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    Z_2_BI_k_1 = structure(c(`1` = 1, `2` = 1, `3` = 1, `4` = 1, 
    `5` = 1, `6` = 1, `7` = 1, `8` = 1, `9` = 1, `10` = 1, `11` = 1, 
    `12` = 1, `13` = 1, `14` = 1, `15` = 1, `16` = 1, `17` = 1, 
    `18` = 1, `19` = 1, `20` = 1, `21` = 1, `22` = 1, `23` = 1, 
    `24` = 1, `25` = 1, `26` = 1, `27` = 1, `28` = 1, `29` = 1, 
    `30` = 1, `31` = 1, `32` = 1, `33` = 1, `34` = 1, `35` = 1
    ), dim = 35L, dimnames = list(c("1", "2", "3", "4", "5", 
    "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", 
    "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
    "27", "28", "29", "30", "31", "32", "33", "34", "35"))), 
    K_BI_int = 1L, X_BI_int = structure(c(1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    K_BI_a = 1L, X_BI_a = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1), dim = c(35L, 1L), dimnames = list(
        c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", 
        "11", "12", "13", "14", "15", "16", "17", "18", "19", 
        "20", "21", "22", "23", "24", "25", "26", "27", "28", 
        "29", "30", "31", "32", "33", "34", "35"), "Intercept"), assign = 0L), 
    J_2_BI = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 
    7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L), dim = 35L), 
    N_1 = 9L, M_1 = 1L, NC_1 = 0L, N_2 = 9L, M_2 = 1L, NC_2 = 0L, 
    prior_only = 0L), class = c("standata", "list"))

And the cod to fit the model:

stan_fit <- stan(file = mod_1,
            data = stan_data,
            warmup = 2000,
            iter = 8000,
            chains = 2,
            cores = 2,
            seed = 17,
            control = list(adapt_delta = 0.9),
            thin = 1)

What can I do to fix this issue? Thank you!

This is very likely due to multimodality. Check the results from different chains and if some chains sample from a infeasible region improve priors. Alternatively, you can use Pathfinder or chain stacking to figure if some of the modes have much less mass.

It was effectively multimodality as I also saw by looking more closely the chains and parameters estimates. I ended-up dropping this analysis as it was not very important for my paper. Thank you for the help!

1 Like