# Parameterization Problems Lead To Poor Mixing

I’m currently modelling the pharmacokinetics of a blood thinner. My goal is to estimate relevant PK parameters through a hierarchical regression.

When using two chains, the chains appear to mix well (all the R-hat statistics are less than 1.1). When I increase the chains to four, the chains stop mixing well, and a lot of parameters see their R-hat increase beyond 2.

The goal of this post is to determine if this is something that I can overcome my drawing more samples, or if there is a fundamental error in my parameterization of the model.

The model can be found here. I’ve tried my best to comment the code well, but feel free to ask more questions as necessary.

Here is the result of `stan_rdump`.

``````D <- 2.5
times <-
c(0.5, 1, 2, 4, 6, 8, 10, 12)
N_t <- 8
N_patients <- 36
C_hat <-
structure(c(7e-04, 0.0082, 0.0133, 0.0156, 0.1021, 0.0337, 0.0374,
0.0263, 0.0823, 0.0066, 0.0752, 0.0067, 0.0043, 0.0106, 0.0046,
0.0159, 0.0179, 0.0039, 0.0071, 0.0282, 0.0211, 0.018, 0.0138,
0.0052, 0.0379, 0.0074, 0.0156, 0.0241, 0.0194, 0.0311, 0.0143,
0.0049, 0.01, 0.0021, 0.0124, 0.0017, 0.0327, 0.0454, 0.1179, 0.0978,
0.1658, 0.0845, 0.0744, 0.0655, 0.1034, 0.0077, 0.0761, 0.0384,
0.0346, 0.0611, 0.0354, 0.0591, 0.0822, 0.0211, 0.0823, 0.0661,
0.0499, 0.0652, 0.041, 0.0432, 0.1199, 0.0183, 0.0392, 0.0505, 0.041,
0.0659, 0.0412, 0.0226, 0.0192, 0.0135, 0.0363, 0.0096, 0.0797,
0.0697, 0.1027, 0.103, 0.1515, 0.1172, 0.0799, 0.162, 0.1041, 0.0129,
0.0771, 0.0944, 0.0838, 0.0678, 0.0322, 0.0781, 0.1265, 0.0428,
0.0787, 0.0641, 0.0468, 0.123, 0.0398, 0.0822, 0.1257, 0.0233,
0.0574, 0.0861, 0.0597, 0.0821, 0.0609, 0.0558, 0.0227, 0.0378,
0.0457, 0.0271, 0.0592, 0.0756, 0.2016, 0.1062, 0.144, 0.0986,
0.0725, 0.1263, 0.081, 0.0448, 0.0735, 0.0838, 0.0728, 0.0671,
0.0271, 0.0758, 0.1055, 0.0661, 0.0585, 0.0616, 0.0425, 0.1053,
0.0421, 0.0803, 0.0984, 0.0297, 0.0686, 0.0914, 0.0713, 0.0903,
0.0583, 0.0418, 0.0323, 0.0796, 0.057, 0.0439, 0.0422, 0.07, 0.119,
0.0774, 0.0942, 0.0652, 0.0536, 0.0936, 0.0524, 0.0378, 0.0646,
0.0793, 0.0447, 0.0417, 0.0167, 0.0431, 0.0652, 0.0525, 0.039,
0.0432, 0.0324, 0.0575, 0.0313, 0.0562, 0.0594, 0.0224, 0.0455,
0.0579, 0.0559, 0.0623, 0.0388, 0.033, 0.03, 0.0741, 0.0428, 0.037,
0.0234, 0.0401, 0.0892, 0.0501, 0.0627, 0.0494, 0.0378, 0.0736,
0.0354, 0.0275, 0.0482, 0.0573, 0.0341, 0.0314, 0.0126, 0.0346,
0.0375, 0.0477, 0.0274, 0.0311, 0.0216, 0.0359, 0.0211, 0.0382,
0.0423, 0.0174, 0.0391, 0.0431, 0.0429, 0.043, 0.0277, 0.0242,
0.0305, 0.061, 0.0394, 0.0278, 0.0164, 0.0268, 0.0659, 0.042, 0.0356,
0.0372, 0.0311, 0.0623, 0.0314, 0.0231, 0.0336, 0.0441, 0.0262,
0.0261, 0.0089, 0.025, 0.0278, 0.0383, 0.0231, 0.029, 0.0174, 0.0225,
0.017, 0.0278, 0.033, 0.0158, 0.0183, 0.033, 0.0355, 0.0306, 0.0199,
0.0208, 0.0254, 0.0507, 0.0376, 0.0218, 0.0329, 0.021, 0.0496,
0.0334, 0.0484, 0.0254, 0.0214, 0.0539, 0.0209, 0.016, 0.0249,
0.0321, 0.0237, 0.0169, 0.0047, 0.0196, 0.0188, 0.0243, 0.0193,
0.0196, 0.0112, 0.0138, 0.0087, 0.018, 0.0225, 0.0108, 0.0123, 0.02,
0.0238, 0.021, 0.0144, 0.0137, 0.0175, 0.0419, 0.0269, 0.0146),
.Dim = c(36, 8))
N_covariates <- 7
X <-
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, 1, 1, 0, 0, 0, 0, 0,
1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0,
0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1,
0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,
-1.00722171115923, -0.563724767446083, -1.30288634030133,
-1.11398949390499, -0.925092647508648, -1.35627023167421,
-0.0339923068998242, -1.13862821300016, -1.15505402573028,
-0.744408707477365, -0.329656936041922, -1.36858959122179,
-0.432318265605151, -0.625321565184021, 1.99459556526957,
-0.157185902375698, 0.0399238503857002, 1.71125029567506,
0.380759464535619, -0.469276344247913, -1.06881850889717,
0.0974141949411084, 0.877640299621645, -0.214676246931107,
-0.797792598850245, 1.4484372919932, 1.96174393980934,
0.298630400885036, 0.877640299621645, 0.816043501883708,
1.37452113470767, 1.05011133328787, 0.820149955066237,
0.766766063693358, 0.701062812772892, -0.411785999692505, -19, 3,
-18, -9, 12, -12, -6, -4, 8, -1, 5, -7, -10, 8, -4, 10, 14, -8, 13,
-18, 20, -15, 0, -24, 15, 1, -11, -8, -1, 0, 13, 7, 10, 17, 2, 11,
1.44244051645102, 2.16480738875956, -0.885186072098717,
-0.403608157226358, -0.162819199790178, -0.483871143038417,
2.16480738875956, -0.243082185602238, -0.243082185602238,
-0.243082185602238, -0.564134128850477, -0.804923086286657,
-0.644397114662537, -0.965449057910776, 0.318758715082181, 0,
-1.44702697278314, -1.44702697278314, 0.72007364414248,
-1.44702697278314, 0.80033662995454, 0.4792846867063,
-0.243082185602238, 1.6832294738872, 1.44244051645102,
-0.483871143038417, 0.72007364414248, -1.28650100115902,
0.63981065833042, 0.80033662995454, -0.323345171414298,
-0.885186072098717, 0.4792846867063, 0.63981065833042,
0.158232743458061, -1.44702697278314, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1,
0, 0, 1),
.Dim = c(36, 7))
A <- 0.03835
``````

## Question

Is my model parameterized poorly, resulting in poor mixing of chains? Can the poor mixing be improved by running a larger warmup?

Your time and help is appreciated.

Hm, nothing sticks out as wrong or parameterized poorly in your model, but I only took a quick look. Since it’s hierarchical you could always be running in to a funnel-like problem, but your priors seems to be pretty informative it seems.

Have your tried looking at pairs plots with the `pairs()` function in Rstan?