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.