Multiple modes on random effect

Modeling
1

Hi,

When trying to fit random effects by looking at the pairs plot I saw several modes. There are bunch of parameters (seven) but most visible is random effect for parameter 2 and 5 (see pair for u.2.4 and u.5.4). This is posterior of parameters 2 & 5 for subject 4. I model random effect for parameter 2 as u2sigma2 where u2~N(0,1) and random effect for parameter 5 as u5sigma5 where u5~N(0,1) . My question is several fold.

Is this approach correct? I “borrowed” it from Matt trick.
Why those modes occur? Is it some nonidentifiability? Or something else?
How to prevent those modes from occurring? Or does it matter?

Can this be caused by assumption that parameters are independent? I would like to get some insight before complicating model with correlation between parameters. It actually exists. The model is below.

Thank you,
Linas

image.png813×815 57.5 KB

2

The simplest answer here is that you’re looking at multiple chains that have not converged. Did you look at the multi-chain R-hat?

3

Also to the general question: yes, even some very simple random effects models are multi-modal.

4

Yes, there is a possibility. Rhat is 1.01594 and neff is 225.573 for u.2.4 and 1.01654 202.948 for u.5.4. I had 4 chains 1200 iterations. Usually Rhat is 1.00… and neff > 1000. Also I see two peaks on the density.

5

Maybe this is the case. Is it bad? Basically I am interested in doing predictions for new individuals. In my model theta_ind=thetaexp(usigma) so new individual will be theta*exp(normal_rng(0,1)*sigma).

6

I would just do the pairs plot and color by chain to see if the chains
converged to different modes

7

If you’re not moving between modes then yes it’s not good in the general
sense. If your chains are not converging you might get sigma wrong as well
so that affects individual predictions.

8

It looks like I am moving between nodes. I think there is nonidentifiability caused by u.2.4*sigma.2. I wonder how to model random effects to remove nonidentifiability.

9

have you looked at traceplots of those to make sure they switch modes often enough? If they’re not it might show up in R-hat or not but a traceplot can tell you qualitatively. It may not be a problem, assuming the model fits, is coded correctly, etc…

10

They change modes pretty often (u.2.4 and u.5.4). But also please take a look at sigma_ind.2 (in the next reply since I am allowed one image per post). My impression that u.2.4 and sigma_ind.2 change modes in sync. Therefore I am thinking that centered parameterization will remove modes. I don’t know if it will affect stan’s performance. I guess it is happening with all u and sigma - only modes are very close.

image.png842×842 205 KB

11

See previous reply.

12

You’re on the right track, when the modes switch this rarely I wouldn’t thing of the chain for sigma_ind.2 as mixing—you might have a rough estimate of the volume of each mode but the n_eff would be low (not sure if that would show up in the estimated n_eff). This might be a situation where the centered parameterization does better since you do have a pretty tight estimate of sigma_ind.2, of course I haven’t read your model so I can’t suggest much more than that.

13

In case you are interested the model is below - I still didn’t figure add how to add code nicely). n_eff is indeed significantly lower than for other parameters. The line with multiple modes is in bold.

functions {
real [] pbpk (real t, real [] y, real [] theta, real [] x_r, int [] x_i) {
real ydot[34];

  real Kpu_scale;
  real Kid_Vmax;
  real Kid_Km;
  real Vmax_ACAT_duo;
  real Vmax_ACAT_jej;
  real Vmax_ACAT_ile;
  real Km_ACAT_Abs;
  real CL_glom;
  real k_leave_1;
  real k_leave_2;
  Kpu_scale = theta[1];
  Kid_Vmax = theta[2];
  Kid_Km = theta[3];
  Vmax_ACAT_duo = theta[4];
  Vmax_ACAT_jej = theta[5];
  Vmax_ACAT_ile = theta[6];
  Km_ACAT_Abs = theta[7];
  CL_glom = theta[8];
  k_leave_1 = theta[9];
  k_leave_2 = theta[10];
  //Vein-1
  //Artery-2
  //Adispose-3
  //Bone-4
  //Brain-5
  //Heart-6
  //Kidney-7
  //Liver-8
  //Lung-9
  //Muscle-10
  //Skin-11
  //Thymus-12
  //Spleen-13
  //Pancreas-14
  //dis-15-23
  //wall-24-32
  //metKidney-33
  //elimColon-34

 ydot[0+1] = 314.4/3.9347*((1./(16.6632+56.9064+16.6632+16.6632+40.2432+13.5192+5.0304+85.5168+63.1944)*(16.6632*y[2+1]*Kpu_scale/(0.335738777114959*0.97/1.1)+56.9064*y[9+1]*Kpu_scale/(1.45080002635366*0.97/1.1)+16.6632*y[10+1]*Kpu_scale/(1.27933750693781*0.97/1.1)+16.6632*y[3+1]*Kpu_scale/(0.753762048673023*0.97/1.1)+40.2432*y[4+1]*Kpu_scale/(0.965738822635123*0.97/1.1)+13.5192*y[5+1]*Kpu_scale/(1.80952944204737*0.97/1.1)+5.0304*y[11+1]*Kpu_scale/(1.83249661975125*0.97/1.1)+85.5168*y[7+1]*Kpu_scale/(2.82861979504438*0.97/1.1)+63.1944*y[6+1]*Kpu_scale/(3.06651237192868*0.97/1.1)))-y[0+1]);
 ydot[1+1] = 314.4/1.9712*((y[8+1]*Kpu_scale/(2.57804385527433*0.97/1.1))-y[1+1]);
 ydot[2+1] = 1./15.4*(16.6632*(y[1+1]-y[2+1]*Kpu_scale/(0.335738777114959*0.97/1.1)));
 ydot[3+1] = 1./11.088*(16.6632*(y[1+1]-y[3+1]*Kpu_scale/(0.753762048673023*0.97/1.1)));
 ydot[4+1] = 1./1.5323*(40.2432*(y[1+1]-y[4+1]*Kpu_scale/(0.965738822635123*0.97/1.1)));
 ydot[5+1] = 1./0.34804*(13.5192*(y[1+1]-y[5+1]*Kpu_scale/(1.80952944204737*0.97/1.1)));
 ydot[6+1] = 1./0.32725*(63.1944*(y[1+1]-y[6+1]*Kpu_scale/(3.06651237192868*0.97/1.1))-(y[6+1]/3.06651237192868*CL_glom+y[6+1]/3.06651237192868*Kid_Vmax/(y[6+1]/3.06651237192868*Kpu_scale+Kid_Km)));
 ydot[7+1] = 1./1.9019*(85.5168*((1./(50.304+10.0608+3.4584+21.6936)*((+3.4584*y[23+1]*Kpu_scale/(1.80953185291572*0.97/1.1)+4.25434530612245*y[24+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[25+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[26+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[27+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[28+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[29+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+4.25434530612245*y[30+1]*Kpu_scale/(1.86402516093248*0.97/1.1)+17.0651828571429*y[31+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+10.0608*y[12+1]*Kpu_scale/(2.2467210773554*0.97/1.1)+3.4584*y[13+1]*Kpu_scale/(1.55110887801799*0.97/1.1)+21.6936*y[1+1]))-y[7+1]*Kpu_scale/(2.82861979504438*0.97/1.1))-0);
 ydot[8+1] = 1./0.52745*(314.4*(y[0+1]-(y[8+1]*Kpu_scale/(2.57804385527433*0.97/1.1))));
 ydot[9+1] = 1./30.569*(56.9064*(y[1+1]-y[9+1]*Kpu_scale/(1.45080002635366*0.97/1.1)));
 ydot[10+1] = 1./3.4804*(16.6632*(y[1+1]-y[10+1]*Kpu_scale/(1.27933750693781*0.97/1.1)));
 ydot[11+1] = 1./0.026334*(5.0304*(y[1+1]-y[11+1]*Kpu_scale/(1.83249661975125*0.97/1.1)));
 ydot[12+1] = 1./0.15785*(10.0608*(y[1+1]-y[12+1]*Kpu_scale/(2.2467210773554*0.97/1.1)));
 ydot[13+1] = 1./0.14784*(3.4584*(y[1+1]-y[13+1]*Kpu_scale/(1.55110887801799*0.97/1.1)));
 ydot[14+1] = 1./0.3*(-k_leave_1*y[14+1]*0.3+0-(0)-0);
 ydot[15+1] = 1./0.105*(k_leave_1*y[14+1]*0.3-k_leave_2*y[15+1]*0.105+0+0-(Vmax_ACAT_duo*y[15+1]/(Km_ACAT_Abs+y[15+1])*0.105)-0);
 ydot[16+1] = 1./0.11*(k_leave_2*y[15+1]*0.105-2.1021021021021*y[16+1]*0.11+0-(Vmax_ACAT_jej*y[16+1]/(Km_ACAT_Abs+y[16+1])*0.11)-0);
 ydot[17+1] = 1./0.11*(2.1021021021021*y[16+1]*0.11-2.1021021021021*y[17+1]*0.11+0-(Vmax_ACAT_jej*y[17+1]/(Km_ACAT_Abs+y[17+1])*0.11)-0);
 ydot[18+1] = 1./0.069*(2.1021021021021*y[17+1]*0.11-2.1021021021021*y[18+1]*0.069+0-(Vmax_ACAT_ile*y[18+1]/(Km_ACAT_Abs+y[18+1])*0.069)-0);
 ydot[19+1] = 1./0.069*(2.1021021021021*y[18+1]*0.069-2.1021021021021*y[19+1]*0.069+0-(Vmax_ACAT_ile*y[19+1]/(Km_ACAT_Abs+y[19+1])*0.069)-0);
 ydot[20+1] = 1./0.069*(2.1021021021021*y[19+1]*0.069-2.1021021021021*y[20+1]*0.069+0-(Vmax_ACAT_ile*y[20+1]/(Km_ACAT_Abs+y[20+1])*0.069)-0);
 ydot[21+1] = 1./0.069*(2.1021021021021*y[20+1]*0.069-2.1021021021021*y[21+1]*0.069+0-(Vmax_ACAT_ile*y[21+1]/(Km_ACAT_Abs+y[21+1])*0.069)-0);
 ydot[22+1] = 1./1*(2.1021021021021*y[21+1]*0.069-0.0840336134453781*y[22+1]*1+0-(0)-0);
 ydot[23+1] = 1./0.147*(3.4584*(y[1+1]-y[23+1]*Kpu_scale/(1.80953185291572*0.97/1.1))+(0)-0-0);
 ydot[24+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[24+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_duo*y[15+1]/(Km_ACAT_Abs+y[15+1])*0.105)-0-0);
 ydot[25+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[25+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_jej*y[16+1]/(Km_ACAT_Abs+y[16+1])*0.11)-0-0);
 ydot[26+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[26+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_jej*y[17+1]/(Km_ACAT_Abs+y[17+1])*0.11)-0-0);
 ydot[27+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[27+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_ile*y[18+1]/(Km_ACAT_Abs+y[18+1])*0.069)-0-0);
 ydot[28+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[28+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_ile*y[19+1]/(Km_ACAT_Abs+y[19+1])*0.069)-0-0);
 ydot[29+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[29+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_ile*y[20+1]/(Km_ACAT_Abs+y[20+1])*0.069)-0-0);
 ydot[30+1] = 1./0.089*(4.25434530612245*(y[1+1]-y[30+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(Vmax_ACAT_ile*y[21+1]/(Km_ACAT_Abs+y[21+1])*0.069)-0-0);
 ydot[31+1] = 1./0.357*(17.0651828571429*(y[1+1]-y[31+1]*Kpu_scale/(1.86402516093248*0.97/1.1))+(0)-0-0);
 ydot[32+1] = (y[6+1]/3.06651237192868*CL_glom+y[6+1]/3.06651237192868*Kid_Vmax/(y[6+1]/3.06651237192868*Kpu_scale+Kid_Km));
 ydot[33+1] = 0.0840336134453781*1*y[22+1];


  return ydot;

}

}

data {

int <lower = 1> T; // number of unique time points

int <lower = 1> N; // number of states

int <lower = 1> Nparam; // number of parameters to estimate

int <lower = 1> Nobs; // number of observations

int <lower = 1> N_IND; // number of individuals

int <lower = 1> N_grp1; // number of individuals in group 1

int <lower = 1, upper = N> idx; // index of state to predict

int <lower = 1, upper = N_IND> ind_idx [Nobs];

// individual for each observation

int <lower = 1, upper = T> ts_idx [Nobs];

// time index for each observation

vector <lower = 0, upper = 1> [N_IND] group_idx;

// group index for each individual

int <lower = 0, upper = N_IND> grp1_idx [N_IND];

// index for individual in group 1

real <lower = 0> y0 [N]; // initial concentrations

real <lower = 0> t0; // initial time

real <lower = 0> uts [T]; // unique time points

vector <lower = 0> [Nobs] y; // observations

vector <lower = 0> [N_IND] urine; // urine observations

int <lower = 1, upper = N> idx_urine; // index of metKidney

}

transformed data {

int x_i [0];

real x_r [0];

vector [N_IND] lurine;

vector [N_IND] one;

vector [Nobs] ly;

int n_raneff;

lurine = log (urine);

one = to_vector (rep_array (1., N_IND));

ly = log (y);

n_raneff = 7;

}

parameters {

real <lower = 0> sigma_prop;

real stheta [Nparam];

real <lower = 0> sigma_urine;

matrix [n_raneff, N_IND] u;

vector [N_grp1] uu;

vector <lower = 0> [n_raneff] sigma_ind;

real <lower = 0> sigma_ind_kid;

}

transformed parameters {

real y_hat_ind [N_IND, T, N];

real <lower = 0> theta_ind [N_IND, Nparam];

vector [Nobs] ly_p;

real <lower = 0> theta [Nparam];

theta [1] = exp (0.1 + 0.06 * stheta [1]); // kpu_scale

theta [2] = exp (5.96 + 2.13 * stheta [2]); // kid_vmax

theta [3] = exp (6.38 + 2.06 * stheta [3]); // kid_km

theta [4] = exp (4.47 + 1.37 * stheta [4]); // acat_vmax_duo

theta [5] = exp (12.48 + 0.94 * stheta [5]); // acat_vmax_jej

theta [6] = exp (7.26 + 1.28 * stheta [6]); // acat_vmax_ile

theta [7] = exp (12.65 + 0.94 * stheta [7]); // acat_km

theta [8] = exp (2.48 + 0.04 * stheta [8]); // CL_glom

theta [9] = exp (0.67 + 0.09 * stheta [9]); // k_leave_1

theta [10] = exp (0.87 + 0.15 * stheta [10]); // k_leave_2

// random effect: log(theta_ind)~N(log(theta),sigma_ind)

theta_ind [, 1] = to_array_1d (theta [1] * exp (u [3] * sigma_ind [3])); // kpu_scale

for (i in 1 : N_IND) {

  if (group_idx [i] == 0) 


     theta_ind [i, 2] = theta [2] * exp (uu [grp1_idx [i]] * sigma_ind_kid); 

     // kid_vmax, group 1


  else


     theta_ind [i, 2] = 0; // kid_vmax=0, group 2

}

theta_ind [, 3] = to_array_1d (theta [3] * one); // kid_km

theta_ind [, 4] = to_array_1d (theta [4] * exp (u [7] * sigma_ind [7])); // acat_vmax_duo

theta_ind [, 5] = to_array_1d (theta [5] * exp (u [2] * sigma_ind [2])); // acat_vmax_jej

theta_ind [, 6] = to_array_1d (theta [6] * exp (u [5] * sigma_ind [5])); // acat_vmax_ile

theta_ind [, 7] = to_array_1d (theta [7] * one); // acat_km

theta_ind [, 8] = to_array_1d (theta [8] * exp (u [4] * sigma_ind [4])); // CL_glom

theta_ind [, 9] = to_array_1d (theta [9] * exp (u [1] * sigma_ind [1])); // k_leave_1

theta_ind [, 10] = to_array_1d (theta [10] * exp (u [6] * sigma_ind [6])); // k_leave_2

for (i in 1 : N_IND)

  y_hat_ind [i] = integrate_ode_bdf (pbpk, y0, t0, uts, theta_ind [i], x_r, x_i);

for (n in 1 : Nobs)

  ly_p [n] = log (y_hat_ind [ind_idx [n], ts_idx [n], idx]);

}

model {

sigma_prop ~ lognormal (-0.99, 0.1);

stheta [1] ~ normal (0, 2); // kpu_scale

stheta [2] ~ normal (0, 2); // kid_vmax

stheta [3] ~ normal (0, 2); // kid_km

stheta [4] ~ normal (0, 2); // acat_vmax_duo

stheta [5] ~ normal (0, 2); // acat_vmax_jej

stheta [6] ~ normal (0, 2); // acat_vmax_ile

stheta [7] ~ normal (0, 2); // acat_km

stheta [8] ~ normal (0, 2); // CL_glom

stheta [9] ~ normal (0, 2); // k_leave_1

stheta [10] ~ normal (0, 2); // k_leave_2

to_vector (u) ~ normal (0, 1);

uu ~ normal (0, 1);

sigma_ind [1] ~ lognormal (-1.03 , 0.7); // k_leave_1

sigma_ind [2] ~ lognormal (-2.72 , 1); // acat_vmax_jej

sigma_ind [3] ~ lognormal (-2.26 , 0.79); // kpu_scale

sigma_ind [4] ~ lognormal (-2.12 , 0.25); // CL_glom

sigma_ind [5] ~ lognormal (-0.02 , 1.02); // acat_vmax_ile

sigma_ind [6] ~ lognormal (-0.86 , 0.7); // k_leave_2

sigma_ind [7] ~ lognormal (0.1 , 1.04); // acat_vmax_duo

sigma_ind_kid ~ lognormal (-0.82 , 1.06); // kid_vmax

ly ~ normal (ly_p, sigma_prop);

// sigma_urine ~ normal (0.26, 0.1); // exp err model

// lurine ~ normal (log (to_vector (y_hat_ind [, T, idx_urine])), sigma_urine); // exp err model

sigma_urine ~ normal (10, 1); // add err model

urine ~ normal (to_vector (y_hat_ind [, T, idx_urine]), sigma_urine); // add err model

}

generated quantities {

real y_rep [N_IND, T];

for (i in 1 : N_IND) {

  real tmp [T, N];



  tmp = integrate_ode_bdf (pbpk, y0, t0, uts, theta_ind [i], x_r, x_i);

  for (t in 1 : T) {

     real lmu;



     lmu = tmp [t, idx];

     if (lmu > 1e-12)

        lmu = log (lmu);

     else

        lmu = -20;

     y_rep [i, t] = exp (normal_rng (lmu, sigma_prop));

  }

}

}