diff --git a/index.html b/index.html
index 17ecd87..417a22f 100644
--- a/index.html
+++ b/index.html
@@ -10,7 +10,7 @@
 <meta name="author" content="Björn Holzhauer - bjoern.holzhauer@novartis.com">
 <meta name="author" content="Lukas Widmer - lukas_andreas.widmer@novartis.com">
 <meta name="author" content="Andrew Bean - andrew.bean@novartis.com">
-<meta name="dcterms.date" content="2024-04-18">
+<meta name="dcterms.date" content="2024-04-19">
 
 <title>(Bayesian) Applied Modelling in Drug Development (BAMDD)</title>
 <style>
@@ -303,7 +303,7 @@ <h2 id="toc-title">Table of contents</h2>
     <div>
     <div class="quarto-title-meta-heading">Published</div>
     <div class="quarto-title-meta-contents">
-      <p class="date">April 18, 2024</p>
+      <p class="date">April 19, 2024</p>
     </div>
   </div>
   
diff --git a/search.json b/search.json
index 6aa7bec..cdd8508 100644
--- a/search.json
+++ b/search.json
@@ -110,7 +110,7 @@
     "href": "src/01b_basic_workflow.html#fitting-the-model",
     "title": "2  Basic workflow",
     "section": "2.5 Fitting the model",
-    "text": "2.5 Fitting the model\nHaving defined the model formula, the family and the priors we can now fit the models with the brm command.\n\nfit_meta_fixed  &lt;- brm(model_meta_fixed, data=arm_data, prior=prior_meta_fixed,\n                       ## setup Stan sampler to be more conservative, but more robust\n                       control=list(adapt_delta=0.95),\n                       ## these options silence Stan\n                       refresh=0, silent=TRUE,\n                       seed=4658758)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.0 seconds.\nChain 2 finished in 0.0 seconds.\nChain 3 finished in 0.0 seconds.\nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.6 seconds.\n\nfit_meta_random  &lt;- brm(model_meta_random, data=arm_data, prior=prior_meta_random,\n                        control=list(adapt_delta=0.95),\n                        refresh=0, silent=TRUE,\n                        seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.9 seconds.\n\n\nWhen calling brm a Stan model file will be created, compiled and run with the data provided. The arguments used are:\n\nformula is the first argument here. It specifies the model. Since we have in this case provided family as part of the formula, we do not need to specify it as argument to brm, which is also possible to do.\ndata the data used to fit.\nprior definition of the priors for all model parameters\ncontrol defines additional control arguments to the Stan sampler. A higher than standard adapt_delta (defaults to 0.8) causes the sampler to run less aggressively which makes the sampler more robust, but somewhat slower.\nrefresh & silent are used here to suppress Stan sampling output.\nseed sets the seed use for the fit.\n\nIn the course of an exploratory analysis one often wishes to modify a given model slightly to study, for example, the sensitivity to different assumptions. For this purpose the update mechanism is a convenient way to simply change certain arguments specified previously. To change the prior on the study level random effect parameter one may use:\n\nprior_meta_random_alt &lt;- prior_meta_fixed +\n    prior(normal(0,0.5), class=\"sd\", coef=\"Intercept\", group=\"study\")\n    \nfit_meta_random_alt  &lt;- update(fit_meta_random, prior=prior_meta_random_alt,\n                               control=list(adapt_delta=0.95),\n                               refresh=0, silent=TRUE,\n                               seed=6845736)\n\nThe desired updates require recompiling the model\n\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.9 seconds.\n\n\nAs changing the prior results in this case in a need to recompile the model, the benefits of using update are not large in this case. However, whenever a model recompilation is not needed, then using update significantly speeds up the workflow as the compilation step is avoided, e.g. when looking at data subsets:\n\nprior_meta_random_alt &lt;- prior_meta_fixed +\n    prior(normal(0,0.5), class=\"sd\", coef=\"Intercept\", group=\"study\")\n    \nfit_meta_random_alt2  &lt;- update(fit_meta_random, newdata=slice_head(arm_data, n=4),\n                                control=list(adapt_delta=0.95),\n                                refresh=0, silent=TRUE,\n                                seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.7 seconds.\n\n\nWarning: 1 of 4000 (0.0%) transitions ended with a divergence.\nSee https://mc-stan.org/misc/warnings for details.\n\n\nNow the fit starts instantly leading to a faster and more convenient modeling workflow.\nWhen working with brms you will encounter warnings now and then which originate from Stan itself. While brms attempts to add explanations to these warnings as to what they mean and how to avoid them, the Stan community has provided an online help-page on possible warnings occuring during a Stan fit.",
+    "text": "2.5 Fitting the model\nHaving defined the model formula, the family and the priors we can now fit the models with the brm command.\n\nfit_meta_fixed  &lt;- brm(model_meta_fixed, data=arm_data, prior=prior_meta_fixed,\n                       ## setup Stan sampler to be more conservative, but more robust\n                       control=list(adapt_delta=0.95),\n                       ## these options silence Stan\n                       refresh=0, silent=TRUE,\n                       seed=4658758)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.0 seconds.\nChain 2 finished in 0.0 seconds.\nChain 3 finished in 0.0 seconds.\nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.6 seconds.\n\nfit_meta_random  &lt;- brm(model_meta_random, data=arm_data, prior=prior_meta_random,\n                        control=list(adapt_delta=0.95),\n                        refresh=0, silent=TRUE,\n                        seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.9 seconds.\n\n\nWhen calling brm a Stan model file will be created, compiled and run with the data provided. The arguments used are:\n\nformula is the first argument here. It specifies the model. Since we have in this case provided family as part of the formula, we do not need to specify it as argument to brm, which is also possible to do.\ndata the data used to fit.\nprior definition of the priors for all model parameters\ncontrol defines additional control arguments to the Stan sampler. A higher than standard adapt_delta (defaults to 0.8) causes the sampler to run less aggressively which makes the sampler more robust, but somewhat slower.\nrefresh & silent are used here to suppress Stan sampling output.\nseed sets the seed use for the fit.\n\nIn the course of an exploratory analysis one often wishes to modify a given model slightly to study, for example, the sensitivity to different assumptions. For this purpose the update mechanism is a convenient way to simply change certain arguments specified previously. To change the prior on the study level random effect parameter one may use:\n\nprior_meta_random_alt &lt;- prior_meta_fixed +\n    prior(normal(0,0.5), class=\"sd\", coef=\"Intercept\", group=\"study\")\n    \nfit_meta_random_alt  &lt;- update(fit_meta_random, prior=prior_meta_random_alt,\n                               control=list(adapt_delta=0.95),\n                               refresh=0, silent=TRUE,\n                               seed=6845736)\n\nThe desired updates require recompiling the model\n\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.9 seconds.\n\n\nAs changing the prior results in this case in a need to recompile the model, the benefits of using update are not large in this case. However, whenever a model recompilation is not needed, then using update significantly speeds up the workflow as the compilation step is avoided, e.g. when looking at data subsets:\n\nprior_meta_random_alt &lt;- prior_meta_fixed +\n    prior(normal(0,0.5), class=\"sd\", coef=\"Intercept\", group=\"study\")\n    \nfit_meta_random_alt2  &lt;- update(fit_meta_random, newdata=slice_head(arm_data, n=4),\n                                control=list(adapt_delta=0.95),\n                                refresh=0, silent=TRUE,\n                                seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.7 seconds.\n\n\nWarning: 1 of 4000 (0.0%) transitions ended with a divergence.\nSee https://mc-stan.org/misc/warnings for details.\n\n\nNow the fit starts instantly leading to a faster and more convenient modeling workflow.\nWhen working with brms you will encounter warnings now and then which originate from Stan itself. While brms attempts to add explanations to these warnings as to what they mean and how to avoid them, the Stan community has provided an online help-page on possible warnings occuring during a Stan fit.",
     "crumbs": [
       "Getting started",
       "<span class='chapter-number'>2</span>  <span class='chapter-title'>Basic workflow</span>"
@@ -121,7 +121,7 @@
     "href": "src/01b_basic_workflow.html#model-summary",
     "title": "2  Basic workflow",
     "section": "2.6 Model summary",
-    "text": "2.6 Model summary\nOnce models are fit, we obtain - by default - a posterior sample of the model. Given that priors are used in such an analysis it can be convenient to first consider the priors which were used to create a given fitting object like:\n\nkable(prior_summary(fit_meta_random))\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\nnormal(0, 2)\nIntercept\n\n\n\n\n\n\n\nuser\n\n\nstudent_t(3, 0, 2.5)\nsd\n\n\n\n\n\n0\n\ndefault\n\n\n\nsd\n\nstudy\n\n\n\n\n\ndefault\n\n\nnormal(0, 1)\nsd\nIntercept\nstudy\n\n\n\n\n\nuser\n\n\n\n\n\nAn even more explicit way to analyze the prior of a model is to sample it directly:\n\nprior_meta_random  &lt;- update(fit_meta_random, sample_prior=\"only\",\n                             control=list(adapt_delta=0.95),\n                             refresh=0, silent=TRUE,\n                             seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.0 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.0 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\n\nThis will sample the model in the usual way, but simply leave out the terms implied by the data likelihood. The resulting model sample can be analyzed in exactly the same way as the model posterior itself. This technique is called prior (predictive) checks.\nReturning to the model posterior, the obvious first thing to do is to simply print the results:\n\nprint(fit_meta_random)\n\n Family: binomial \n  Links: mu = logit \nFormula: r | trials(n) ~ 1 + (1 | study) \n   Data: arm_data (Number of observations: 8) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nMultilevel Hyperparameters:\n~study (Number of levels: 8) \n              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsd(Intercept)     0.39      0.21     0.04     0.89 1.00     1040      969\n\nRegression Coefficients:\n          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept    -1.10      0.19    -1.46    -0.71 1.00     1491     1834\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nWe see that we fitted the default 4 chains and ran each chain for 1000 warmup iterations and 2000 total iterations such that we are left with 4k iterations overall. brms reports by default as estimate the mean of the posterior with it’s standard error and the 95% credible interval. The Rhat column is a convergence diagnostic which should be close to 1.0. Values above 1.1 indicate non-convergence of the Markov Chain for the particular parameter. The additional two columns on bulk and tail ESS indicate statistical information on central moments and tail properties of the target quantitiy. The ESS is the effective sample size of the MC estimator, which assumes that the estimation error of the mean scales with \\(1/\\sqrt{ESS}\\). An ESS of ~200 is usually sufficient for a precise estimate of the mean. It is important to note that the Stan HMC sampler often requires fewer iterations (for which more time is spent) as compared to BUGS or JAGS in order to reach a comparable ESS. For more details on the ESS, please refer to the Stan reference manual on ESS.\nbrms supports the common R functions to extract model results:\n\nfitted to get the posterior estimates of expected values of for each data row (no sampling uncertainty)\npredict to get the posterior predictive estimates of the response for each data row (includes sampling uncertainty)\ncoef/ranef to get the random effect estimates including/excluding the linear predictor part\nmany more, please refer to the reference manual of brms\n\nOf note is the argument summary, which is used in a number of these functions. The default is set to TRUE such that the output are summaries of the posterior sample in the form of means, credible intervals, etc. When setting the argument summary=FALSE then the posterior sample is returned in the format of a matrix. Each row of the matrix is one iteration while the columns of the matrix run with the rows of the input data-set.\nFor example, we can compare the estimated mean response fo the two different models as:\n\nfitted(fit_meta_random)\n\n      Estimate Est.Error      Q2.5     Q97.5\n[1,] 24.331054  3.735007 17.183654 31.875063\n[2,] 11.544028  2.191310  7.565441 16.378010\n[3,] 16.002830  2.961917 11.039539 22.567141\n[4,]  9.412544  1.869250  5.841321 13.179245\n[5,] 37.769373  4.748738 29.170767 47.811069\n[6,]  5.359016  1.232845  3.174211  8.080594\n[7,] 13.531353  3.402142  7.355087 20.241300\n[8,]  9.335581  1.922931  6.012867 13.677123\n\nfitted(fit_meta_fixed)\n\n      Estimate Est.Error      Q2.5     Q97.5\n[1,] 26.533687 2.0380505 22.574910 30.640824\n[2,] 10.911049 0.8380768  9.283141 12.599965\n[3,] 12.646898 0.9714072 10.760004 14.604505\n[4,]  9.671157 0.7428408  8.228238 11.168151\n[5,] 34.468995 2.6475609 29.326285 39.804435\n[6,]  4.959568 0.3809440  4.219609  5.727257\n[7,] 19.342314 1.4856817 16.456477 22.336302\n[8,]  8.679243 0.6666520  7.384316 10.022699\n\n\nAs expected, the standard errors for the fixed effect analysis are considerably smaller.",
+    "text": "2.6 Model summary\nOnce models are fit, we obtain - by default - a posterior sample of the model. Given that priors are used in such an analysis it can be convenient to first consider the priors which were used to create a given fitting object like:\n\nkable(prior_summary(fit_meta_random))\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\nnormal(0, 2)\nIntercept\n\n\n\n\n\n\n\nuser\n\n\nstudent_t(3, 0, 2.5)\nsd\n\n\n\n\n\n0\n\ndefault\n\n\n\nsd\n\nstudy\n\n\n\n\n\ndefault\n\n\nnormal(0, 1)\nsd\nIntercept\nstudy\n\n\n\n\n\nuser\n\n\n\n\n\nAn even more explicit way to analyze the prior of a model is to sample it directly:\n\nprior_meta_random  &lt;- update(fit_meta_random, sample_prior=\"only\",\n                             control=list(adapt_delta=0.95),\n                             refresh=0, silent=TRUE,\n                             seed=5868467)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.0 seconds.\nChain 2 finished in 0.0 seconds.\nChain 3 finished in 0.0 seconds.\nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\n\nThis will sample the model in the usual way, but simply leave out the terms implied by the data likelihood. The resulting model sample can be analyzed in exactly the same way as the model posterior itself. This technique is called prior (predictive) checks.\nReturning to the model posterior, the obvious first thing to do is to simply print the results:\n\nprint(fit_meta_random)\n\n Family: binomial \n  Links: mu = logit \nFormula: r | trials(n) ~ 1 + (1 | study) \n   Data: arm_data (Number of observations: 8) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nMultilevel Hyperparameters:\n~study (Number of levels: 8) \n              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsd(Intercept)     0.39      0.21     0.04     0.89 1.00     1040      969\n\nRegression Coefficients:\n          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept    -1.10      0.19    -1.46    -0.71 1.00     1491     1834\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nWe see that we fitted the default 4 chains and ran each chain for 1000 warmup iterations and 2000 total iterations such that we are left with 4k iterations overall. brms reports by default as estimate the mean of the posterior with it’s standard error and the 95% credible interval. The Rhat column is a convergence diagnostic which should be close to 1.0. Values above 1.1 indicate non-convergence of the Markov Chain for the particular parameter. The additional two columns on bulk and tail ESS indicate statistical information on central moments and tail properties of the target quantitiy. The ESS is the effective sample size of the MC estimator, which assumes that the estimation error of the mean scales with \\(1/\\sqrt{ESS}\\). An ESS of ~200 is usually sufficient for a precise estimate of the mean. It is important to note that the Stan HMC sampler often requires fewer iterations (for which more time is spent) as compared to BUGS or JAGS in order to reach a comparable ESS. For more details on the ESS, please refer to the Stan reference manual on ESS.\nbrms supports the common R functions to extract model results:\n\nfitted to get the posterior estimates of expected values of for each data row (no sampling uncertainty)\npredict to get the posterior predictive estimates of the response for each data row (includes sampling uncertainty)\ncoef/ranef to get the random effect estimates including/excluding the linear predictor part\nmany more, please refer to the reference manual of brms\n\nOf note is the argument summary, which is used in a number of these functions. The default is set to TRUE such that the output are summaries of the posterior sample in the form of means, credible intervals, etc. When setting the argument summary=FALSE then the posterior sample is returned in the format of a matrix. Each row of the matrix is one iteration while the columns of the matrix run with the rows of the input data-set.\nFor example, we can compare the estimated mean response fo the two different models as:\n\nfitted(fit_meta_random)\n\n      Estimate Est.Error      Q2.5     Q97.5\n[1,] 24.331054  3.735007 17.183654 31.875063\n[2,] 11.544028  2.191310  7.565441 16.378010\n[3,] 16.002830  2.961917 11.039539 22.567141\n[4,]  9.412544  1.869250  5.841321 13.179245\n[5,] 37.769373  4.748738 29.170767 47.811069\n[6,]  5.359016  1.232845  3.174211  8.080594\n[7,] 13.531353  3.402142  7.355087 20.241300\n[8,]  9.335581  1.922931  6.012867 13.677123\n\nfitted(fit_meta_fixed)\n\n      Estimate Est.Error      Q2.5     Q97.5\n[1,] 26.533687 2.0380505 22.574910 30.640824\n[2,] 10.911049 0.8380768  9.283141 12.599965\n[3,] 12.646898 0.9714072 10.760004 14.604505\n[4,]  9.671157 0.7428408  8.228238 11.168151\n[5,] 34.468995 2.6475609 29.326285 39.804435\n[6,]  4.959568 0.3809440  4.219609  5.727257\n[7,] 19.342314 1.4856817 16.456477 22.336302\n[8,]  8.679243 0.6666520  7.384316 10.022699\n\n\nAs expected, the standard errors for the fixed effect analysis are considerably smaller.",
     "crumbs": [
       "Getting started",
       "<span class='chapter-number'>2</span>  <span class='chapter-title'>Basic workflow</span>"
@@ -373,7 +373,7 @@
     "href": "src/02b_dose_finding.html#implementation-and-results",
     "title": "5  Dose finding",
     "section": "5.6 Implementation and results",
-    "text": "5.6 Implementation and results\nWe will now fit the various dose response models introduced in the statistical model section using brms. We will start with the sigmoid-Emax model and look at the things brms let us do easily with a single fitted model. We will not repeat that level of exploration for other models, but you may wish to do so in practice.\n\n5.6.1 Fitting dose response models with brms: sigmoid-Emax model\nWe use the capabilities of brms for fitting non-linear models. For example, a sigmoid-Emax model can be fitted as follows:\n\n# Fitting a sigmoid Emax model\nbrmfit1 &lt;-  brm( bf( est | se(stderr) ~ E0 + Emax * dose^h/(dose^h + ED50^h),\n                    nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),\n                    E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, \n                    nl=TRUE, family=gaussian()),\n                data=pathway,\n                prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                          prior(normal(0,1), nlpar=\"logh\"),\n                          prior(normal(0,1), nlpar=\"Emax\"),\n                          prior(normal(4,2), nlpar=\"logED50\")),\n                control=list(adapt_delta=0.999),\n                seed=3624)\n\nNote that with individual patient data, we would not write est | se(stderr) ~, but would instead use outcome ~ and bf() + negbinomial() (and similarly, e.g. + gaussian() for continuous data). With the settings above, you might still experience a divergent transition or two and you may wish to increase adapt_delta to 0.9999 or 0.99999, or alternatively try to come up with a way to reparameterize the model.\n\n# Summarize model fit\nsummary(brmfit1)\n\nWarning: There were 1 divergent transitions after warmup. Increasing\nadapt_delta above 0.999 may help. See\nhttp://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup\n\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * dose^h/(dose^h + ED50^h) \n         h ~ exp(logh)\n         ED50 ~ exp(logED50)\n         E0 ~ 1\n         logED50 ~ 1\n         logh ~ 1\n         Emax ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept         -0.42      0.10    -0.61    -0.21 1.00     2050     2044\nlogED50_Intercept     2.73      1.29    -0.10     4.99 1.00     1251     1246\nlogh_Intercept        0.00      0.95    -1.82     1.85 1.00     1397     1605\nEmax_Intercept       -1.28      0.30    -2.05    -0.83 1.00     1502     1168\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSo what does the model fit look like? We have two options for looking at this. The easiest and fastest uses the conditional_effects function from brms.\n\nplot(conditional_effects(brmfit1, \"dose\"), \n     points = TRUE)\n\n\n\n\n\n\n\n\nThe downside is that in a controlled trial we are primarily interested in the contrast with placebo and that we would like to see the confidence intervals around the point estimates for each dose. We can obtain that with a little bit more code:\n\n# We are going to predict the annualized rate for all doses from 0,1,2,...600.\n# We are setting stderr=0 to avoid predicting with sampling variability.\ndr_curve_data = tibble(dose=as.integer(seq(0,600)), stderr=0)\n\n# Now we get the predictions and turn this into a dataset with one row per dose per MCMC sample\ndr_curve_pred &lt;- tidybayes::add_predicted_draws(dr_curve_data, brmfit1) %&gt;%\n  ungroup() %&gt;%\n  transmute(dose, .sample = .draw, pred = .prediction)\n\n# Now we get the rate ratio by subtracting the annualized placebo rate for each MCMC sample  \ndr_curve_pred = dr_curve_pred %&gt;%\n  left_join(dr_curve_pred  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;%\n  dplyr::select(dose, .sample, logRR)\n\n# Now we get median prediction, 50 and 95% credible intervals \n# and plot those with the original data overlaid.\np1 = dr_curve_pred  %&gt;%\n  dplyr::select(-.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"red\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"red\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkred\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"SigMod Emax model fit\"), col=\"darkred\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(p1)\n\n\n\n\n\n\n\n\nWe can also plot the dose response curve for each MCMC sample:\n\n dr_curve_pred  %&gt;%\n  ggplot(aes(x=dose, y=logRR, group=factor(.sample))) +\n  geom_hline(yintercept=0)  +\n  geom_line(col=\"darkred\", alpha=0.1, size=0.25) +\n  geom_point(data=pathway_deltas %&gt;% mutate(.sample=NA_integer_)) +\n  geom_errorbar(data=pathway_deltas %&gt;% mutate(.sample=NA_integer_), \n                aes(ymin=.lower, ymax=.upper), width=10) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\n\n\n\n\n\n\n\nNote, that this also makes it really easy to determine a distribution for the lowest dose that achieves a 50% relative rate reduction:\n\ndr_curve_pred  %&gt;%\n  filter(logRR&lt;log(0.5) | (logRR&gt;=log(0.5) & dose==600)) %&gt;%\n  group_by(.sample) %&gt;%\n  summarize(dose = min(dose), logRR=max(logRR)) %&gt;%\n  mutate(`Has 50% RRR`=(logRR&lt;log(0.5))) %&gt;%\n  ggplot(aes(x=dose, y=\"Dose with\\n&gt;50% RRR\")) +\n  stat_halfeye() +\n  coord_cartesian(xlim=c(0,200)) +\n  xlab(\"Tezepelumab dose [mg/month]\") +\n  ylab(\"\")\n\n\n\n\n\n\n\n\nAs the plot below shows the data resulted in a substantially different posterior distribution compared with the prior distribution for the logED50 parameter and the Emax parameter. These parameters are reasonably well-identified even with just 3 doses - although the ED50 would be better estimated, if doses with approximately 50% of the Emax (maximum relative rate reduction) had been included in the study. In contrast, the log-Hill-parameter is poorly identified, when there are very few doses (as in this study), which is why the posterior distribution is more or less the prior distribution for the logED50 parameter.\n\nbrmfit1 %&gt;%\n  spread_draws(b_logED50_Intercept, b_logh_Intercept, b_Emax_Intercept) %&gt;%\n  pivot_longer(cols=c(\"b_logED50_Intercept\", \"b_logh_Intercept\", \"b_Emax_Intercept\")) %&gt;%\n  ggplot(aes(x=value)) +\n  theme_bw(base_size=16) +\n  geom_density() +\n  geom_line(data=tibble(name=c(rep(\"b_logED50_Intercept\",151),\n                               rep(\"b_logh_Intercept\",71),\n                               rep(\"b_Emax_Intercept\",61)),\n                        xval = c(seq(-5,10,0.1),\n                                 seq(-3,4,0.1),\n                                 seq(-3,3,0.1)),\n                        density=\n                          case_when(str_detect(name,\"ED50\") ~ dnorm(xval,mean=4,sd=2),\n                                    str_detect(name,\"logh\") ~  dnorm(xval,mean=0,sd=1),\n                                    TRUE ~ dnorm(xval,mean=0,sd=1))),\n            aes(x=xval, y=density), col=\"red\") +\n  geom_text(data=tibble(name=\"b_Emax_Intercept\", density=0.6, xval=1),\n            aes(x=xval, y=density, label=\"prior\"), col=\"red\", size=8) +\n  geom_text(data=tibble(name=\"b_Emax_Intercept\", density=1, xval=0.5),\n            aes(x=xval, y=density, label=\"posterior\"), col=\"black\", size=8) +\n  facet_wrap(~name, scales = \"free\")\n\n\n\n\n\n\n\n\n\n\n5.6.2 Fitting dose response models with brms: modified beta model\nWe now fit the modified beta model that can capture non-monotone dose response relationships.\n\n# Try a modified beta model\nbrmfit2 &lt;-  brm( bf( est | se(stderr) ~ E0 + Emax*(delta1+delta2)^(delta1+delta2)/(delta1^delta1*delta2^delta2)*(dose/850)^delta1*(1-dose/850)^delta2,\n                   E0 ~ 1, delta1 ~ 1, delta2 ~ 1, Emax ~ 1,\n                   nl=TRUE, \n                   family = gaussian()),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(lognormal(0,1), nlpar=\"delta1\", lb=0),\n                         prior(lognormal(0,1), nlpar=\"delta2\", lb=0)),\n               control=list(adapt_delta=0.999),\n               save_pars = save_pars(all = TRUE),\n               seed = 7304)\n\nLet us summarize the model fit:\n\nsummary(brmfit2)\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * (delta1 + delta2)^(delta1 + delta2)/(delta1^delta1 * delta2^delta2) * (dose/850)^delta1 * (1 - dose/850)^delta2 \n         E0 ~ 1\n         delta1 ~ 1\n         delta2 ~ 1\n         Emax ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept        -0.43      0.10    -0.63    -0.23 1.00     1535     2122\ndelta1_Intercept     0.45      0.22     0.13     0.98 1.00     1303     1812\ndelta2_Intercept     0.75      0.42     0.16     1.75 1.00     1304     1577\nEmax_Intercept      -1.30      0.21    -1.73    -0.90 1.00     1640     1769\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSo what does the model fit look like? We could, again, just use plot(conditional_effects(brmfit2, \"dose\"), points = TRUE) or get a slightly more tailored plot.\n\n# Now we get the predictions and turn this into a dataset with one row per dose per MCMC sample\ndr_curve_pred2 &lt;- tidybayes::add_predicted_draws(dr_curve_data, brmfit2) %&gt;%\n  ungroup() %&gt;%\n  transmute(dose, .sample = .draw, pred = .prediction)\n\n# Now we get the rate ratio by subtracting the annualized placebo rate for each MCMC sample  \ndr_curve_pred2 = dr_curve_pred2 %&gt;%\n  left_join(dr_curve_pred2  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;%\n  dplyr::select(dose, .sample, logRR)\n\n# Now we get median prediction, 50 and 95% credible intervals \n# and plot those with the original data overlaid.\np2 = dr_curve_pred2  %&gt;%\n  dplyr::select(-.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"royalblue\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"royalblue\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkblue\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"modified-beta model fit\"), col=\"darkblue\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(p2)\n\n\n\n\n\n\n\n\n\n\n5.6.3 Perform Bayesian model averaging\nNow we have two fitted model. We can see both of them together below:\n\np1 + p2\n\n\n\n\n\n\n\n\nIn practice, we do not really know, whether one of these is the true data generating model. In fact, realistically both models are going to be at least somewhat misspecified, because we made no attempt to fully model the full complexity of the underlying biological system. So, really, the question is which of the two models provides the better approximation. At a glance both seem to fit the data somewhat decently, so it seems unlikely that we can completely rule one of the two models out. That is where Bayesian model averaging comes in.\nFirst, we need to compare how well each model fits. As it turns out, when you just have 4 data points (due to us using summary data), it is problematic to just use the “standard” brms::loo function to perform approximate leave-one-out cross-validation (LOO-CV) based on the posterior likelihood as implemented in the loo package, as can be seen by running\n\nloo(brmfit1)\n\nWarning: Found 1 observations with a pareto_k &gt; 0.7 in model 'brmfit1'. We\nrecommend to set 'moment_match = TRUE' in order to perform moment matching for\nproblematic observations.\n\n\n\nComputed from 4000 by 4 log-likelihood matrix.\n\n         Estimate  SE\nelpd_loo      1.0 0.3\np_loo         1.3 0.5\nlooic        -2.1 0.7\n------\nMCSE of elpd_loo is NA.\nMCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.9]).\n\nPareto k diagnostic values:\n                         Count Pct.    Min. ESS\n(-Inf, 0.7]   (good)     3     75.0%   852     \n   (0.7, 1]   (bad)      1     25.0%   &lt;NA&gt;    \n   (1, Inf)   (very bad) 0      0.0%   &lt;NA&gt;    \nSee help('pareto-k-diagnostic') for details.\n\n\nand\n\nloo(brmfit2)\n\nWarning: Found 4 observations with a pareto_k &gt; 0.7 in model 'brmfit2'. We\nrecommend to set 'moment_match = TRUE' in order to perform moment matching for\nproblematic observations.\n\n\n\nComputed from 4000 by 4 log-likelihood matrix.\n\n         Estimate  SE\nelpd_loo     -0.4 0.4\np_loo         2.7 0.7\nlooic         0.7 0.8\n------\nMCSE of elpd_loo is NA.\nMCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.6]).\n\nPareto k diagnostic values:\n                         Count Pct.    Min. ESS\n(-Inf, 0.7]   (good)     0      0.0%   &lt;NA&gt;    \n   (0.7, 1]   (bad)      2     50.0%   &lt;NA&gt;    \n   (1, Inf)   (very bad) 2     50.0%   &lt;NA&gt;    \nSee help('pareto-k-diagnostic') for details.\n\n\nInstead, we actually fit each model leaving each dose group out. We could probably do something more sophisticated by simulating consistent individual patient data (and then performing approximate LOO-CV) or parametric sampling from the normal distribution around the point estimates.\n\n(loo_exact_brmfit1 &lt;- kfold(brmfit1, folds = \"loo\"))\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.6 seconds.\nChain 2 finished in 0.8 seconds.\nChain 3 finished in 0.8 seconds.\nChain 4 finished in 0.8 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.7 seconds.\nTotal execution time: 3.3 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.2 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.5 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.8 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.4 seconds.\nChain 3 finished in 0.4 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.8 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.7 seconds.\n\n\n\nBased on 4-fold cross-validation.\n\n           Estimate  SE\nelpd_kfold     -0.5 1.3\np_kfold         2.8 1.9\nkfoldic         1.0 2.7\n\n(loo_exact_brmfit2 &lt;- kfold(brmfit2, folds = \"loo\"))\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.6 seconds.\nChain 2 finished in 0.5 seconds.\nChain 3 finished in 0.4 seconds.\nChain 4 finished in 0.4 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.5 seconds.\nTotal execution time: 2.3 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.2 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.5 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.4 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.6 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.5 seconds.\nChain 3 finished in 0.5 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.4 seconds.\nTotal execution time: 2.0 seconds.\n\n\n\nBased on 4-fold cross-validation.\n\n           Estimate  SE\nelpd_kfold     -2.1 1.7\np_kfold         4.4 2.2\nkfoldic         4.2 3.4\n\n\nNow, we can compare the models via\n\nloo_compare(loo_exact_brmfit1, loo_exact_brmfit2)\n\n        elpd_diff se_diff\nbrmfit1  0.0       0.0   \nbrmfit2 -1.6       0.8   \n\n\n\n5.6.3.1 Approach 1 (manual without brms)\nLet us say that we a-priori assign a 75% probability to the SigEmax model and 25% to the modified beta model.\n\n# follows: \n# Gould, A. L. (2019). BMA‐Mod: A Bayesian model averaging strategy for determining dose‐response relationships in the presence of model uncertainty. Biometrical Journal, 61(5), 1141-1159.\n# https://dx.doi.org/10.1002/bimj.201700211\nprior_weights = c(0.75,0.25)\nposterior_weigths = c(prior_weights[1]*exp(loo_exact_brmfit1$estimates[\"elpd_kfold\", \"Estimate\"]), \n                      prior_weights[2]*exp(loo_exact_brmfit2$estimates[\"elpd_kfold\", \"Estimate\"]))\n(posterior_weigths = posterior_weigths/sum(posterior_weigths))\n\n[1] 0.93629621 0.06370379\n\n\nAs we can see, using leaving-one-dose-group out cross-validation, suggests that the SigEmax model should be given most of the weight. A-posteriori, we assign 93.6% weight to the SigEmax model and 6.4% weight to the modified beta model. However, if we had done approximate LOO-CV with individual patient data, the modified beta model might have been given more weight.\nWhen we average the predictions of the two models with these weights, we get the dose response curve below.\n\nbma_curve = dr_curve_pred %&gt;%\n  inner_join(dr_curve_pred2, by=c(\".sample\", \"dose\")) %&gt;%\n  mutate(logRR = posterior_weigths[1]*logRR.x + posterior_weigths[2]*logRR.y) %&gt;%\n  dplyr::select(-logRR.x, -logRR.y, -.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95))\n\n\npbma = bma_curve %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"orange\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"orange\", alpha=0.5) +\n  geom_line(col=\"darkorange\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"Bayesian model averaging\"), col=\"darkorange\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose\")\n\nplot( ( (p1 + theme(axis.title.x=element_blank()) + ggtitle(\"Sigmoid Emax model\")) | \n          (p2 + theme(axis.title.y=element_blank(),\n                      axis.title.x=element_blank()) + \n             ggtitle(\"Modified beta model\"))) / \n        (pbma + ggtitle(\"Bayesian model averaging\") | plot_spacer()) )\n\n\n\n\n\n\n\n\n\n\n5.6.3.2 Approach 2: Using brms functions\nBefore, we needed to do a lot of manual work to get model-averaged predictions, but the below using brms::pp_average would be a lot simpler, if we could rely on the approximate LOO-CV.\n\n5.6.3.2.1 Approach 2a: Not appropriate here\n\nbma_brms_curve = dr_curve_data %&gt;%\n  dplyr::select(-stderr) %&gt;%\n  bind_cols(as_tibble(\n    t(pp_average(brmfit1, brmfit2, newdata=dr_curve_data, \n                 weights=\"loo\", # Not a good option in this case.\n                 summary=F)),\n    .name_repair=vctrs::vec_as_names(repair=\"unique\", quiet =T))) %&gt;%\n  pivot_longer(cols=starts_with(\"V\"), names_to=\".sample\", values_to=\"pred\") %&gt;%\n  mutate(.sample = as.integer(str_extract(.sample, \"[0-9]+\")))\n\n# Forming differences to placebo by merging with the placebo (dose=0) predictions for each sample\nbma_brms_curve = bma_brms_curve %&gt;%\n  left_join(bma_brms_curve  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;% # Calculate log rate ratio (difference on log scale vs. placebo)\n  dplyr::select(dose, .sample, logRR)\n\n\n\n5.6.3.2.2 Approach 2b: Using leave-one-dose-out CV\nBecause we want to avoid using approximate LOO-CV in this case, we first replace the internally calculated approximate LOO-CV results with our own\n\nbrmfit1$criteria$loo &lt;- loo_exact_brmfit1\nbrmfit2$criteria$loo &lt;- loo_exact_brmfit2\n(w_dose &lt;- model_weights(brmfit1, brmfit2, weights = \"loo\"))\n\n brmfit1  brmfit2 \n0.830486 0.169514 \n\n\nWe can then do the following:\n\nbmapreds1 &lt;- posterior_epred(brmfit1, newdata = dr_curve_data)\nbmapreds2 &lt;- posterior_epred(brmfit2, newdata = dr_curve_data)\npe_avg &lt;- bmapreds1 * w_dose[1] + bmapreds2 * w_dose[2]\nrownames(pe_avg) &lt;- 1:nrow(pe_avg)\n\n# log-rate predictions for each dose\npe_avg = pe_avg %&gt;% \n  t() %&gt;% \n  as_tibble() %&gt;%\n  bind_cols(dr_curve_data %&gt;% dplyr::select(dose)) %&gt;%\n  pivot_longer(cols=-dose, names_to = \".sample\", values_to=\"pred\") \n\n# form difference to placebo (log-rate-ratio) and plot\n pbma2 = pe_avg %&gt;%\n  left_join(y=pe_avg %&gt;% filter(dose==0) %&gt;% dplyr::select(-dose), \n            by=\".sample\") %&gt;%\n  mutate(logRR=pred.x-pred.y) %&gt;%\n  dplyr::select(-.sample, -pred.x, -pred.y) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"green\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"darkgreen\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkgreen\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"BMA (LOO weights)\"), col=\"darkgreen\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(pbma2)\n\n\n\n\n\n\n\n\n\n\n\n5.6.3.3 Why did we do this form of Bayesian model averaging?\nWhy are we doing it like this? Can we not just fit a single model like this (as an example where we average a sigmoid-Emax and a modified beta model)? \\[\n\\begin{aligned}\nf(\\text{dose}; \\text{parameters}) = & \\text{E}_0 + w \\times \\text{E}_{\\text{max }1} \\times \\frac{\\text{dose}^h}{\\text{dose}^h + \\text{ED}_{50}^h} \\\\\n  & + (1-w) \\times \\text{E}_{\\text{max }2} \\times \\frac{(\\delta_1 + \\delta_2)^{\\delta_1 + \\delta_2}}{\\delta_1^{\\delta_1} + \\delta_2^{\\delta_2}} \\times (\\frac{\\text{dose}}{S})^{\\delta_1} \\times (1-\\frac{\\text{dose}}{S})^{\\delta_2}\n\\end{aligned}\n\\] brms certainly lets us specify such a model (see below). But, part of the reason why we do not use it, is that the model parameters are really badly identified: you can still obtain a similarly good fit at each observed dose level by making one of the two models fit better for one dose (and worse for another dose), if the other model in turn is made to fit worse for that dose (and better for the other dose).\nIf you run the code below, you will see how much the sampler struggles to fit this model, if we try to fit it. That is the case, even if we make things easier and fix the model weights to 50:50.\n\n# The model below does not work so well\nbrmfit3 &lt;-  brm( bf( est | se(stderr) ~ E0 + 0.5 * Emax * dose^h/(dose^h + ED50^h) + 0.5*Emax2*(delta1+delta2)^(delta1+delta2)/(delta1^delta1*delta2^delta2)*(dose/850)^delta1*(1-dose/850)^delta2,\n                   E0 ~ 1, ED50 ~ 1, h ~ 1, Emax ~ 1,\n                   #modelwgt ~ 1,\n                   delta1 ~ 1, delta2 ~ 1, Emax2 ~ 1,\n                   nl=T),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(lognormal(0,1), nlpar=\"h\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(lognormal(4,2), nlpar=\"ED50\"),\n                         #prior(beta(4,1), nlpar=\"modelwgt\", lb=0, ub=1),\n                         prior(normal(0,1), nlpar=\"Emax2\"),\n                         prior(lognormal(0,1), nlpar=\"delta1\", lb=0),\n                         prior(lognormal(0,1), nlpar=\"delta2\", lb=0)),\n               control=list(adapt_delta=0.999))\n\nIf you run this code, it produces a lot of severe warnings (which you should not ignore) about divergent transitions after warmup and a very high largest R-hat. Additionally, “Bulk Effective Samples Size” and “Tail Effective Samples Size (ESS)” are too low, and a lot of transitions after warmup that exceeded the maximum treedepth.\n\n\n\n5.6.4 Going beyond the default MCP-Mod models\nSo far, a lot of what we did could also be done in a frequentist manner using the MCP-Mod approach. However, with brms it is easy to make our models more complex when necessary. For example, we assumed that a 280 mg q2w dose is equivalent to a 560 mg q4w dose. This is of course only the case in terms of the total amount of drug injected over the trial period. However, the two dose regimens differ in terms of their blood concentrations over time as shown in the plot below.\n\nsim_1st_order_pk = function(injected=c(1,rep(0,28*5)), dose=1, thalf=28*24, V=10){\n  k = log(2)/thalf  \n  concentration = rep(0, length(injected))\n  for (i in 1:length(injected)){\n    if (i&gt;1) {\n      concentration[i] = concentration[i-1] - k*concentration[i-1] + dose*injected[i]/V\n    } else {\n      concentration[i] = dose*injected[i]/V\n    }\n  }\n  return(concentration)\n}\n\nplot_data = tibble(regimen=\"q4wk\",\n                   dose=280,\n                   concentration = sim_1st_order_pk(rep(c(1,rep(0, 28*24-1)),12), dose=280)) %&gt;%\n  mutate(hour=1:n()) %&gt;%\n  bind_rows(tibble(regimen=\"q2wk\",\n       dose=280,\n       concentration = sim_1st_order_pk(rep(c(1,rep(0, 14*24-1)),24), dose=280)) %&gt;%\n         mutate(hour=1:n())) %&gt;%\n  bind_rows(bind_rows(tibble(regimen=\"q4wk\",\n       dose=560,\n       concentration = sim_1st_order_pk(rep(c(1,rep(0, 28*24-1)),12), dose=560)) %&gt;%\n         mutate(hour=1:n()))) %&gt;%\n  mutate(dose_regimen = paste0(dose, \" mg \", regimen)) \n\nplot_data %&gt;%\n  ggplot(aes(x=hour, y=concentration, col=dose_regimen)) +\n  geom_line() +\n  coord_cartesian(ylim=c(0,120)) +\n  scale_x_continuous(breaks=seq(0,12)*28*24, labels=seq(0,12)) +\n  geom_text(data=. %&gt;% filter( (hour==8025 & str_detect(dose_regimen, \"280 mg q2\")) |\n                               (hour==500 & str_detect(dose_regimen, \"560 mg q4\")) |\n                               (hour==5000 & str_detect(dose_regimen, \"280 mg q4\"))),\n            aes(label=dose_regimen), nudge_y=c(-20,25, 60), angle=c(0,-75,0)) +\n  xlab(\"Time in trial (months)\") +\n  ylab(\"Drug concentration [mg/L]\")\n\n\n\n\n\n\n\n\nThese differences in pharmacokinetic profiles likely lead to at least somewhat different efficacy outcomes. For example, if efficacy is more driven by the peak concentration achieved in each dosing interval, 280 mg q2w might be less effective than a 560 mg q4w regiment would have been. On the other hand, if efficacy is more determined by the minimum concentration maintained throughout time, then 280 mg q2w could be more effective than 560 mg q4w. If the main determinant of efficacy is the average drug concentration over time, and peaks and troughs do not matter (at least within the concentrations seen), then the two dosing regimens might be exactly identical.\nSo, we could update our previous sigmoid Emax model and add an extra term that describes how the 280 mg q2w dose might relate to q4w doses. Note that assuming a monotone concentration response, a 280 mg q2w regimen cannot be worse than a 280 mg q4w regimen - this is also illustrated by the plot showing that the concentrations of 280 mg q2w are (unsurprisingly) always above those of 280 mg q4w. This gives us a lower bound for the efficacy of 280 mg q2w. Similarly, its peak concentrations stay below 2.5 times the peak concentrations of 280 mg q2w, which gives us an upper bound. We will specify a log-normal prior with mean 0 and SD 0.67 for a regimen factor that indicates what q4w dose a q2w dose is equivalent to.\n\nbrmfit3 &lt;- brm( bf( est | se(stderr) ~ E0 + Emax * (dose*dosefactor)^h / ((dose*dosefactor)^h + ED50^h),\n                   nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)), nlf(dosefactor ~ exp((dose==560)*regimen)),\n                   E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, regimen ~ 1, \n                   nl=T, family=gaussian()),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(normal(0,1), nlpar=\"logh\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(normal(4,2), nlpar=\"logED50\"),\n                         prior(normal(0,0.67), nlpar=\"regimen\")),\n               control=list(adapt_delta=0.999),\n               seed=2023,\n               save_pars = save_pars(all = TRUE))\n\nLet’s summarize the model results:\n\nsummary(brmfit3)\n\nWarning: There were 3 divergent transitions after warmup. Increasing\nadapt_delta above 0.999 may help. See\nhttp://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup\n\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * (dose * dosefactor)^h/((dose * dosefactor)^h + ED50^h) \n         h ~ exp(logh)\n         ED50 ~ exp(logED50)\n         dosefactor ~ exp((dose == 560) * regimen)\n         E0 ~ 1\n         logED50 ~ 1\n         logh ~ 1\n         Emax ~ 1\n         regimen ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept         -0.42      0.10    -0.61    -0.22 1.00     2624     2529\nlogED50_Intercept     2.76      1.34    -0.33     5.09 1.00     1552     1553\nlogh_Intercept        0.00      0.98    -1.81     1.98 1.00     1388     1769\nEmax_Intercept       -1.30      0.33    -2.11    -0.84 1.00     1429     1348\nregimen_Intercept    -0.05      0.68    -1.36     1.25 1.00     2553     2452\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nAs we can see, we can fit such a model, but the marginal posterior for the log-regimen factor is still a N(0, 0.67) distribution. Thus, the data do not inform this parameter very well - which should not surprise us. However, if we had information external to the PATHWAY trial, such as the elicited judgments of experts, we could incorporate it via this type of model.",
+    "text": "5.6 Implementation and results\nWe will now fit the various dose response models introduced in the statistical model section using brms. We will start with the sigmoid-Emax model and look at the things brms let us do easily with a single fitted model. We will not repeat that level of exploration for other models, but you may wish to do so in practice.\n\n5.6.1 Fitting dose response models with brms: sigmoid-Emax model\nWe use the capabilities of brms for fitting non-linear models. For example, a sigmoid-Emax model can be fitted as follows:\n\n# Fitting a sigmoid Emax model\nbrmfit1 &lt;-  brm( bf( est | se(stderr) ~ E0 + Emax * dose^h/(dose^h + ED50^h),\n                    nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),\n                    E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, \n                    nl=TRUE, family=gaussian()),\n                data=pathway,\n                prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                          prior(normal(0,1), nlpar=\"logh\"),\n                          prior(normal(0,1), nlpar=\"Emax\"),\n                          prior(normal(4,2), nlpar=\"logED50\")),\n                control=list(adapt_delta=0.999),\n                seed=3624)\n\nNote that with individual patient data, we would not write est | se(stderr) ~, but would instead use outcome ~ and bf() + negbinomial() (and similarly, e.g. + gaussian() for continuous data). With the settings above, you might still experience a divergent transition or two and you may wish to increase adapt_delta to 0.9999 or 0.99999, or alternatively try to come up with a way to reparameterize the model.\n\n# Summarize model fit\nsummary(brmfit1)\n\nWarning: There were 1 divergent transitions after warmup. Increasing\nadapt_delta above 0.999 may help. See\nhttp://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup\n\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * dose^h/(dose^h + ED50^h) \n         h ~ exp(logh)\n         ED50 ~ exp(logED50)\n         E0 ~ 1\n         logED50 ~ 1\n         logh ~ 1\n         Emax ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept         -0.42      0.10    -0.61    -0.21 1.00     2050     2044\nlogED50_Intercept     2.73      1.29    -0.10     4.99 1.00     1251     1246\nlogh_Intercept        0.00      0.95    -1.82     1.85 1.00     1397     1605\nEmax_Intercept       -1.28      0.30    -2.05    -0.83 1.00     1502     1168\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSo what does the model fit look like? We have two options for looking at this. The easiest and fastest uses the conditional_effects function from brms.\n\nplot(conditional_effects(brmfit1, \"dose\"), \n     points = TRUE)\n\n\n\n\n\n\n\n\nThe downside is that in a controlled trial we are primarily interested in the contrast with placebo and that we would like to see the confidence intervals around the point estimates for each dose. We can obtain that with a little bit more code:\n\n# We are going to predict the annualized rate for all doses from 0,1,2,...600.\n# We are setting stderr=0 to avoid predicting with sampling variability.\ndr_curve_data = tibble(dose=as.integer(seq(0,600)), stderr=0)\n\n# Now we get the predictions and turn this into a dataset with one row per dose per MCMC sample\ndr_curve_pred &lt;- tidybayes::add_predicted_draws(dr_curve_data, brmfit1) %&gt;%\n  ungroup() %&gt;%\n  transmute(dose, .sample = .draw, pred = .prediction)\n\n# Now we get the rate ratio by subtracting the annualized placebo rate for each MCMC sample  \ndr_curve_pred = dr_curve_pred %&gt;%\n  left_join(dr_curve_pred  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;%\n  dplyr::select(dose, .sample, logRR)\n\n# Now we get median prediction, 50 and 95% credible intervals \n# and plot those with the original data overlaid.\np1 = dr_curve_pred  %&gt;%\n  dplyr::select(-.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"red\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"red\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkred\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"SigMod Emax model fit\"), col=\"darkred\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(p1)\n\n\n\n\n\n\n\n\nWe can also plot the dose response curve for each MCMC sample:\n\n dr_curve_pred  %&gt;%\n  ggplot(aes(x=dose, y=logRR, group=factor(.sample))) +\n  geom_hline(yintercept=0)  +\n  geom_line(col=\"darkred\", alpha=0.1, size=0.25) +\n  geom_point(data=pathway_deltas %&gt;% mutate(.sample=NA_integer_)) +\n  geom_errorbar(data=pathway_deltas %&gt;% mutate(.sample=NA_integer_), \n                aes(ymin=.lower, ymax=.upper), width=10) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\n\n\n\n\n\n\n\nNote, that this also makes it really easy to determine a distribution for the lowest dose that achieves a 50% relative rate reduction:\n\ndr_curve_pred  %&gt;%\n  filter(logRR&lt;log(0.5) | (logRR&gt;=log(0.5) & dose==600)) %&gt;%\n  group_by(.sample) %&gt;%\n  summarize(dose = min(dose), logRR=max(logRR)) %&gt;%\n  mutate(`Has 50% RRR`=(logRR&lt;log(0.5))) %&gt;%\n  ggplot(aes(x=dose, y=\"Dose with\\n&gt;50% RRR\")) +\n  stat_halfeye() +\n  coord_cartesian(xlim=c(0,200)) +\n  xlab(\"Tezepelumab dose [mg/month]\") +\n  ylab(\"\")\n\n\n\n\n\n\n\n\nAs the plot below shows the data resulted in a substantially different posterior distribution compared with the prior distribution for the logED50 parameter and the Emax parameter. These parameters are reasonably well-identified even with just 3 doses - although the ED50 would be better estimated, if doses with approximately 50% of the Emax (maximum relative rate reduction) had been included in the study. In contrast, the log-Hill-parameter is poorly identified, when there are very few doses (as in this study), which is why the posterior distribution is more or less the prior distribution for the logED50 parameter.\n\nbrmfit1 %&gt;%\n  spread_draws(b_logED50_Intercept, b_logh_Intercept, b_Emax_Intercept) %&gt;%\n  pivot_longer(cols=c(\"b_logED50_Intercept\", \"b_logh_Intercept\", \"b_Emax_Intercept\")) %&gt;%\n  ggplot(aes(x=value)) +\n  theme_bw(base_size=16) +\n  geom_density() +\n  geom_line(data=tibble(name=c(rep(\"b_logED50_Intercept\",151),\n                               rep(\"b_logh_Intercept\",71),\n                               rep(\"b_Emax_Intercept\",61)),\n                        xval = c(seq(-5,10,0.1),\n                                 seq(-3,4,0.1),\n                                 seq(-3,3,0.1)),\n                        density=\n                          case_when(str_detect(name,\"ED50\") ~ dnorm(xval,mean=4,sd=2),\n                                    str_detect(name,\"logh\") ~  dnorm(xval,mean=0,sd=1),\n                                    TRUE ~ dnorm(xval,mean=0,sd=1))),\n            aes(x=xval, y=density), col=\"red\") +\n  geom_text(data=tibble(name=\"b_Emax_Intercept\", density=0.6, xval=1),\n            aes(x=xval, y=density, label=\"prior\"), col=\"red\", size=8) +\n  geom_text(data=tibble(name=\"b_Emax_Intercept\", density=1, xval=0.5),\n            aes(x=xval, y=density, label=\"posterior\"), col=\"black\", size=8) +\n  facet_wrap(~name, scales = \"free\")\n\n\n\n\n\n\n\n\n\n\n5.6.2 Fitting dose response models with brms: modified beta model\nWe now fit the modified beta model that can capture non-monotone dose response relationships.\n\n# Try a modified beta model\nbrmfit2 &lt;-  brm( bf( est | se(stderr) ~ E0 + Emax*(delta1+delta2)^(delta1+delta2)/(delta1^delta1*delta2^delta2)*(dose/850)^delta1*(1-dose/850)^delta2,\n                   E0 ~ 1, delta1 ~ 1, delta2 ~ 1, Emax ~ 1,\n                   nl=TRUE, \n                   family = gaussian()),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(lognormal(0,1), nlpar=\"delta1\", lb=0),\n                         prior(lognormal(0,1), nlpar=\"delta2\", lb=0)),\n               control=list(adapt_delta=0.999),\n               save_pars = save_pars(all = TRUE),\n               seed = 7304)\n\nLet us summarize the model fit:\n\nsummary(brmfit2)\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * (delta1 + delta2)^(delta1 + delta2)/(delta1^delta1 * delta2^delta2) * (dose/850)^delta1 * (1 - dose/850)^delta2 \n         E0 ~ 1\n         delta1 ~ 1\n         delta2 ~ 1\n         Emax ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept        -0.43      0.10    -0.63    -0.23 1.00     1535     2122\ndelta1_Intercept     0.45      0.22     0.13     0.98 1.00     1303     1812\ndelta2_Intercept     0.75      0.42     0.16     1.75 1.00     1304     1577\nEmax_Intercept      -1.30      0.21    -1.73    -0.90 1.00     1640     1769\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSo what does the model fit look like? We could, again, just use plot(conditional_effects(brmfit2, \"dose\"), points = TRUE) or get a slightly more tailored plot.\n\n# Now we get the predictions and turn this into a dataset with one row per dose per MCMC sample\ndr_curve_pred2 &lt;- tidybayes::add_predicted_draws(dr_curve_data, brmfit2) %&gt;%\n  ungroup() %&gt;%\n  transmute(dose, .sample = .draw, pred = .prediction)\n\n# Now we get the rate ratio by subtracting the annualized placebo rate for each MCMC sample  \ndr_curve_pred2 = dr_curve_pred2 %&gt;%\n  left_join(dr_curve_pred2  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;%\n  dplyr::select(dose, .sample, logRR)\n\n# Now we get median prediction, 50 and 95% credible intervals \n# and plot those with the original data overlaid.\np2 = dr_curve_pred2  %&gt;%\n  dplyr::select(-.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"royalblue\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"royalblue\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkblue\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"modified-beta model fit\"), col=\"darkblue\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(p2)\n\n\n\n\n\n\n\n\n\n\n5.6.3 Perform Bayesian model averaging\nNow we have two fitted model. We can see both of them together below:\n\np1 + p2\n\n\n\n\n\n\n\n\nIn practice, we do not really know, whether one of these is the true data generating model. In fact, realistically both models are going to be at least somewhat misspecified, because we made no attempt to fully model the full complexity of the underlying biological system. So, really, the question is which of the two models provides the better approximation. At a glance both seem to fit the data somewhat decently, so it seems unlikely that we can completely rule one of the two models out. That is where Bayesian model averaging comes in.\nFirst, we need to compare how well each model fits. As it turns out, when you just have 4 data points (due to us using summary data), it is problematic to just use the “standard” brms::loo function to perform approximate leave-one-out cross-validation (LOO-CV) based on the posterior likelihood as implemented in the loo package, as can be seen by running\n\nloo(brmfit1)\n\nWarning: Found 1 observations with a pareto_k &gt; 0.7 in model 'brmfit1'. We\nrecommend to set 'moment_match = TRUE' in order to perform moment matching for\nproblematic observations.\n\n\n\nComputed from 4000 by 4 log-likelihood matrix.\n\n         Estimate  SE\nelpd_loo      1.0 0.3\np_loo         1.3 0.5\nlooic        -2.1 0.7\n------\nMCSE of elpd_loo is NA.\nMCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.9]).\n\nPareto k diagnostic values:\n                         Count Pct.    Min. ESS\n(-Inf, 0.7]   (good)     3     75.0%   852     \n   (0.7, 1]   (bad)      1     25.0%   &lt;NA&gt;    \n   (1, Inf)   (very bad) 0      0.0%   &lt;NA&gt;    \nSee help('pareto-k-diagnostic') for details.\n\n\nand\n\nloo(brmfit2)\n\nWarning: Found 4 observations with a pareto_k &gt; 0.7 in model 'brmfit2'. We\nrecommend to set 'moment_match = TRUE' in order to perform moment matching for\nproblematic observations.\n\n\n\nComputed from 4000 by 4 log-likelihood matrix.\n\n         Estimate  SE\nelpd_loo     -0.4 0.4\np_loo         2.7 0.7\nlooic         0.7 0.8\n------\nMCSE of elpd_loo is NA.\nMCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.6]).\n\nPareto k diagnostic values:\n                         Count Pct.    Min. ESS\n(-Inf, 0.7]   (good)     0      0.0%   &lt;NA&gt;    \n   (0.7, 1]   (bad)      2     50.0%   &lt;NA&gt;    \n   (1, Inf)   (very bad) 2     50.0%   &lt;NA&gt;    \nSee help('pareto-k-diagnostic') for details.\n\n\nInstead, we actually fit each model leaving each dose group out. We could probably do something more sophisticated by simulating consistent individual patient data (and then performing approximate LOO-CV) or parametric sampling from the normal distribution around the point estimates.\n\n(loo_exact_brmfit1 &lt;- kfold(brmfit1, folds = \"loo\"))\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.6 seconds.\nChain 2 finished in 0.7 seconds.\nChain 3 finished in 0.8 seconds.\nChain 4 finished in 0.8 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.7 seconds.\nTotal execution time: 3.2 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.2 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.4 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.5 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.4 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.6 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.7 seconds.\n\n\n\nBased on 4-fold cross-validation.\n\n           Estimate  SE\nelpd_kfold     -0.5 1.3\np_kfold         2.8 1.9\nkfoldic         1.0 2.7\n\n(loo_exact_brmfit2 &lt;- kfold(brmfit2, folds = \"loo\"))\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.6 seconds.\nChain 2 finished in 0.5 seconds.\nChain 3 finished in 0.4 seconds.\nChain 4 finished in 0.4 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.5 seconds.\nTotal execution time: 2.3 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.3 seconds.\nChain 4 finished in 0.2 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.5 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.3 seconds.\nChain 2 finished in 0.3 seconds.\nChain 3 finished in 0.4 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.3 seconds.\nTotal execution time: 1.6 seconds.\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.2 seconds.\nChain 2 finished in 0.5 seconds.\nChain 3 finished in 0.5 seconds.\nChain 4 finished in 0.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.4 seconds.\nTotal execution time: 2.0 seconds.\n\n\n\nBased on 4-fold cross-validation.\n\n           Estimate  SE\nelpd_kfold     -2.1 1.7\np_kfold         4.4 2.2\nkfoldic         4.2 3.4\n\n\nNow, we can compare the models via\n\nloo_compare(loo_exact_brmfit1, loo_exact_brmfit2)\n\n        elpd_diff se_diff\nbrmfit1  0.0       0.0   \nbrmfit2 -1.6       0.8   \n\n\n\n5.6.3.1 Approach 1 (manual without brms)\nLet us say that we a-priori assign a 75% probability to the SigEmax model and 25% to the modified beta model.\n\n# follows: \n# Gould, A. L. (2019). BMA‐Mod: A Bayesian model averaging strategy for determining dose‐response relationships in the presence of model uncertainty. Biometrical Journal, 61(5), 1141-1159.\n# https://dx.doi.org/10.1002/bimj.201700211\nprior_weights = c(0.75,0.25)\nposterior_weigths = c(prior_weights[1]*exp(loo_exact_brmfit1$estimates[\"elpd_kfold\", \"Estimate\"]), \n                      prior_weights[2]*exp(loo_exact_brmfit2$estimates[\"elpd_kfold\", \"Estimate\"]))\n(posterior_weigths = posterior_weigths/sum(posterior_weigths))\n\n[1] 0.93629621 0.06370379\n\n\nAs we can see, using leaving-one-dose-group out cross-validation, suggests that the SigEmax model should be given most of the weight. A-posteriori, we assign 93.6% weight to the SigEmax model and 6.4% weight to the modified beta model. However, if we had done approximate LOO-CV with individual patient data, the modified beta model might have been given more weight.\nWhen we average the predictions of the two models with these weights, we get the dose response curve below.\n\nbma_curve = dr_curve_pred %&gt;%\n  inner_join(dr_curve_pred2, by=c(\".sample\", \"dose\")) %&gt;%\n  mutate(logRR = posterior_weigths[1]*logRR.x + posterior_weigths[2]*logRR.y) %&gt;%\n  dplyr::select(-logRR.x, -logRR.y, -.sample) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95))\n\n\npbma = bma_curve %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"orange\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"orange\", alpha=0.5) +\n  geom_line(col=\"darkorange\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"Bayesian model averaging\"), col=\"darkorange\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose\")\n\nplot( ( (p1 + theme(axis.title.x=element_blank()) + ggtitle(\"Sigmoid Emax model\")) | \n          (p2 + theme(axis.title.y=element_blank(),\n                      axis.title.x=element_blank()) + \n             ggtitle(\"Modified beta model\"))) / \n        (pbma + ggtitle(\"Bayesian model averaging\") | plot_spacer()) )\n\n\n\n\n\n\n\n\n\n\n5.6.3.2 Approach 2: Using brms functions\nBefore, we needed to do a lot of manual work to get model-averaged predictions, but the below using brms::pp_average would be a lot simpler, if we could rely on the approximate LOO-CV.\n\n5.6.3.2.1 Approach 2a: Not appropriate here\n\nbma_brms_curve = dr_curve_data %&gt;%\n  dplyr::select(-stderr) %&gt;%\n  bind_cols(as_tibble(\n    t(pp_average(brmfit1, brmfit2, newdata=dr_curve_data, \n                 weights=\"loo\", # Not a good option in this case.\n                 summary=F)),\n    .name_repair=vctrs::vec_as_names(repair=\"unique\", quiet =T))) %&gt;%\n  pivot_longer(cols=starts_with(\"V\"), names_to=\".sample\", values_to=\"pred\") %&gt;%\n  mutate(.sample = as.integer(str_extract(.sample, \"[0-9]+\")))\n\n# Forming differences to placebo by merging with the placebo (dose=0) predictions for each sample\nbma_brms_curve = bma_brms_curve %&gt;%\n  left_join(bma_brms_curve  %&gt;% \n              filter(dose==0) %&gt;% \n              dplyr::select(-dose) %&gt;% \n              rename(placebo=pred), \n            by=\".sample\") %&gt;%\n  mutate(logRR = pred - placebo) %&gt;% # Calculate log rate ratio (difference on log scale vs. placebo)\n  dplyr::select(dose, .sample, logRR)\n\n\n\n5.6.3.2.2 Approach 2b: Using leave-one-dose-out CV\nBecause we want to avoid using approximate LOO-CV in this case, we first replace the internally calculated approximate LOO-CV results with our own\n\nbrmfit1$criteria$loo &lt;- loo_exact_brmfit1\nbrmfit2$criteria$loo &lt;- loo_exact_brmfit2\n(w_dose &lt;- model_weights(brmfit1, brmfit2, weights = \"loo\"))\n\n brmfit1  brmfit2 \n0.830486 0.169514 \n\n\nWe can then do the following:\n\nbmapreds1 &lt;- posterior_epred(brmfit1, newdata = dr_curve_data)\nbmapreds2 &lt;- posterior_epred(brmfit2, newdata = dr_curve_data)\npe_avg &lt;- bmapreds1 * w_dose[1] + bmapreds2 * w_dose[2]\nrownames(pe_avg) &lt;- 1:nrow(pe_avg)\n\n# log-rate predictions for each dose\npe_avg = pe_avg %&gt;% \n  t() %&gt;% \n  as_tibble() %&gt;%\n  bind_cols(dr_curve_data %&gt;% dplyr::select(dose)) %&gt;%\n  pivot_longer(cols=-dose, names_to = \".sample\", values_to=\"pred\") \n\n# form difference to placebo (log-rate-ratio) and plot\n pbma2 = pe_avg %&gt;%\n  left_join(y=pe_avg %&gt;% filter(dose==0) %&gt;% dplyr::select(-dose), \n            by=\".sample\") %&gt;%\n  mutate(logRR=pred.x-pred.y) %&gt;%\n  dplyr::select(-.sample, -pred.x, -pred.y) %&gt;%\n  group_by(dose) %&gt;%\n  median_qi(.width = c(0.5, 0.95)) %&gt;%\n  ungroup() %&gt;%\n  ggplot(aes(x=dose, y=logRR, ymin=.lower, ymax=.upper)) +\n  geom_hline(yintercept=0)  +\n  geom_ribbon(data=.%&gt;% filter(.width==0.95), fill=\"green\", alpha=0.5) +\n  geom_ribbon(data=.%&gt;% filter(.width==0.5), fill=\"darkgreen\", alpha=0.5) +\n  geom_line(data=.%&gt;% filter(.width==0.5), col=\"darkgreen\") +\n  geom_point(data=pathway_deltas) +\n  geom_errorbar(data=pathway_deltas, width=10) +\n  geom_text(data=tibble(dose=300, logRR=log(0.6), .lower=NA_real_, .upper=NA_real_),\n            aes(label=\"BMA (LOO weights)\"), col=\"darkgreen\", size=4) +\n  scale_y_continuous(breaks=log(c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \n                     labels=paste(100*(1-c(0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1)), \"%\")) +\n  ylab(\"Relative rate reduction\") +\n  xlab(\"Total tezepelumab dose [mg/month]\")\n\nplot(pbma2)\n\n\n\n\n\n\n\n\n\n\n\n5.6.3.3 Why did we do this form of Bayesian model averaging?\nWhy are we doing it like this? Can we not just fit a single model like this (as an example where we average a sigmoid-Emax and a modified beta model)? \\[\n\\begin{aligned}\nf(\\text{dose}; \\text{parameters}) = & \\text{E}_0 + w \\times \\text{E}_{\\text{max }1} \\times \\frac{\\text{dose}^h}{\\text{dose}^h + \\text{ED}_{50}^h} \\\\\n  & + (1-w) \\times \\text{E}_{\\text{max }2} \\times \\frac{(\\delta_1 + \\delta_2)^{\\delta_1 + \\delta_2}}{\\delta_1^{\\delta_1} + \\delta_2^{\\delta_2}} \\times (\\frac{\\text{dose}}{S})^{\\delta_1} \\times (1-\\frac{\\text{dose}}{S})^{\\delta_2}\n\\end{aligned}\n\\] brms certainly lets us specify such a model (see below). But, part of the reason why we do not use it, is that the model parameters are really badly identified: you can still obtain a similarly good fit at each observed dose level by making one of the two models fit better for one dose (and worse for another dose), if the other model in turn is made to fit worse for that dose (and better for the other dose).\nIf you run the code below, you will see how much the sampler struggles to fit this model, if we try to fit it. That is the case, even if we make things easier and fix the model weights to 50:50.\n\n# The model below does not work so well\nbrmfit3 &lt;-  brm( bf( est | se(stderr) ~ E0 + 0.5 * Emax * dose^h/(dose^h + ED50^h) + 0.5*Emax2*(delta1+delta2)^(delta1+delta2)/(delta1^delta1*delta2^delta2)*(dose/850)^delta1*(1-dose/850)^delta2,\n                   E0 ~ 1, ED50 ~ 1, h ~ 1, Emax ~ 1,\n                   #modelwgt ~ 1,\n                   delta1 ~ 1, delta2 ~ 1, Emax2 ~ 1,\n                   nl=T),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(lognormal(0,1), nlpar=\"h\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(lognormal(4,2), nlpar=\"ED50\"),\n                         #prior(beta(4,1), nlpar=\"modelwgt\", lb=0, ub=1),\n                         prior(normal(0,1), nlpar=\"Emax2\"),\n                         prior(lognormal(0,1), nlpar=\"delta1\", lb=0),\n                         prior(lognormal(0,1), nlpar=\"delta2\", lb=0)),\n               control=list(adapt_delta=0.999))\n\nIf you run this code, it produces a lot of severe warnings (which you should not ignore) about divergent transitions after warmup and a very high largest R-hat. Additionally, “Bulk Effective Samples Size” and “Tail Effective Samples Size (ESS)” are too low, and a lot of transitions after warmup that exceeded the maximum treedepth.\n\n\n\n5.6.4 Going beyond the default MCP-Mod models\nSo far, a lot of what we did could also be done in a frequentist manner using the MCP-Mod approach. However, with brms it is easy to make our models more complex when necessary. For example, we assumed that a 280 mg q2w dose is equivalent to a 560 mg q4w dose. This is of course only the case in terms of the total amount of drug injected over the trial period. However, the two dose regimens differ in terms of their blood concentrations over time as shown in the plot below.\n\nsim_1st_order_pk = function(injected=c(1,rep(0,28*5)), dose=1, thalf=28*24, V=10){\n  k = log(2)/thalf  \n  concentration = rep(0, length(injected))\n  for (i in 1:length(injected)){\n    if (i&gt;1) {\n      concentration[i] = concentration[i-1] - k*concentration[i-1] + dose*injected[i]/V\n    } else {\n      concentration[i] = dose*injected[i]/V\n    }\n  }\n  return(concentration)\n}\n\nplot_data = tibble(regimen=\"q4wk\",\n                   dose=280,\n                   concentration = sim_1st_order_pk(rep(c(1,rep(0, 28*24-1)),12), dose=280)) %&gt;%\n  mutate(hour=1:n()) %&gt;%\n  bind_rows(tibble(regimen=\"q2wk\",\n       dose=280,\n       concentration = sim_1st_order_pk(rep(c(1,rep(0, 14*24-1)),24), dose=280)) %&gt;%\n         mutate(hour=1:n())) %&gt;%\n  bind_rows(bind_rows(tibble(regimen=\"q4wk\",\n       dose=560,\n       concentration = sim_1st_order_pk(rep(c(1,rep(0, 28*24-1)),12), dose=560)) %&gt;%\n         mutate(hour=1:n()))) %&gt;%\n  mutate(dose_regimen = paste0(dose, \" mg \", regimen)) \n\nplot_data %&gt;%\n  ggplot(aes(x=hour, y=concentration, col=dose_regimen)) +\n  geom_line() +\n  coord_cartesian(ylim=c(0,120)) +\n  scale_x_continuous(breaks=seq(0,12)*28*24, labels=seq(0,12)) +\n  geom_text(data=. %&gt;% filter( (hour==8025 & str_detect(dose_regimen, \"280 mg q2\")) |\n                               (hour==500 & str_detect(dose_regimen, \"560 mg q4\")) |\n                               (hour==5000 & str_detect(dose_regimen, \"280 mg q4\"))),\n            aes(label=dose_regimen), nudge_y=c(-20,25, 60), angle=c(0,-75,0)) +\n  xlab(\"Time in trial (months)\") +\n  ylab(\"Drug concentration [mg/L]\")\n\n\n\n\n\n\n\n\nThese differences in pharmacokinetic profiles likely lead to at least somewhat different efficacy outcomes. For example, if efficacy is more driven by the peak concentration achieved in each dosing interval, 280 mg q2w might be less effective than a 560 mg q4w regiment would have been. On the other hand, if efficacy is more determined by the minimum concentration maintained throughout time, then 280 mg q2w could be more effective than 560 mg q4w. If the main determinant of efficacy is the average drug concentration over time, and peaks and troughs do not matter (at least within the concentrations seen), then the two dosing regimens might be exactly identical.\nSo, we could update our previous sigmoid Emax model and add an extra term that describes how the 280 mg q2w dose might relate to q4w doses. Note that assuming a monotone concentration response, a 280 mg q2w regimen cannot be worse than a 280 mg q4w regimen - this is also illustrated by the plot showing that the concentrations of 280 mg q2w are (unsurprisingly) always above those of 280 mg q4w. This gives us a lower bound for the efficacy of 280 mg q2w. Similarly, its peak concentrations stay below 2.5 times the peak concentrations of 280 mg q2w, which gives us an upper bound. We will specify a log-normal prior with mean 0 and SD 0.67 for a regimen factor that indicates what q4w dose a q2w dose is equivalent to.\n\nbrmfit3 &lt;- brm( bf( est | se(stderr) ~ E0 + Emax * (dose*dosefactor)^h / ((dose*dosefactor)^h + ED50^h),\n                   nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)), nlf(dosefactor ~ exp((dose==560)*regimen)),\n                   E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, regimen ~ 1, \n                   nl=T, family=gaussian()),\n               data=pathway,\n               prior = c(prior(normal(0,1), nlpar=\"E0\"),\n                         prior(normal(0,1), nlpar=\"logh\"),\n                         prior(normal(0,1), nlpar=\"Emax\"),\n                         prior(normal(4,2), nlpar=\"logED50\"),\n                         prior(normal(0,0.67), nlpar=\"regimen\")),\n               control=list(adapt_delta=0.999),\n               seed=2023,\n               save_pars = save_pars(all = TRUE))\n\nLet’s summarize the model results:\n\nsummary(brmfit3)\n\nWarning: There were 3 divergent transitions after warmup. Increasing\nadapt_delta above 0.999 may help. See\nhttp://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup\n\n\n Family: gaussian \n  Links: mu = identity; sigma = identity \nFormula: est | se(stderr) ~ E0 + Emax * (dose * dosefactor)^h/((dose * dosefactor)^h + ED50^h) \n         h ~ exp(logh)\n         ED50 ~ exp(logED50)\n         dosefactor ~ exp((dose == 560) * regimen)\n         E0 ~ 1\n         logED50 ~ 1\n         logh ~ 1\n         Emax ~ 1\n         regimen ~ 1\n   Data: pathway (Number of observations: 4) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nE0_Intercept         -0.42      0.10    -0.61    -0.22 1.00     2624     2529\nlogED50_Intercept     2.76      1.34    -0.33     5.09 1.00     1552     1553\nlogh_Intercept        0.00      0.98    -1.81     1.98 1.00     1388     1769\nEmax_Intercept       -1.30      0.33    -2.11    -0.84 1.00     1429     1348\nregimen_Intercept    -0.05      0.68    -1.36     1.25 1.00     2553     2452\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nsigma     0.00      0.00     0.00     0.00   NA       NA       NA\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nAs we can see, we can fit such a model, but the marginal posterior for the log-regimen factor is still a N(0, 0.67) distribution. Thus, the data do not inform this parameter very well - which should not surprise us. However, if we had information external to the PATHWAY trial, such as the elicited judgments of experts, we could incorporate it via this type of model.",
     "crumbs": [
       "Case studies",
       "<span class='chapter-number'>5</span>  <span class='chapter-title'>Dose finding</span>"
@@ -527,7 +527,7 @@
     "href": "src/02e_multiple_imputation.html#implementation",
     "title": "7  Multiple imputation",
     "section": "7.5 Implementation",
-    "text": "7.5 Implementation\n\n7.5.1 Imputation from observed off-treatment data\nFirst, let us conduct an imputation model based on the observed off-treatment data, in which we specify that event rates can be different on- and off-treatment for each treatment group. Using this model, we perform predictions for the missing post-treatment values (patients lost-to-follow-up).\nWe use family=poisson() here, because a Poisson model with a random patient effect on the intercept is approximately equivalent to a negative binomial regression model and this is a way of reflecting the correlation of multiple observation periods for a patient that is consistent with that. However, we could also make the imputation model more complex e.g. by allowing counts to be overdispersed within patient by using family=negbinomial(), or allowing for some interactions such as treatment:logBASE or period:logBASE.\n\n# Fit imputation model\nbrmfit1 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + period + treatment:period + logBASE,\n               data=count_data %&gt;% \n                   filter(period %in% c(\"On-treatment\", \"Post-treatment\") & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 11.3 seconds.\nChain 2 finished in 10.6 seconds.\nChain 3 finished in 13.7 seconds.\nChain 4 finished in 13.2 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 12.2 seconds.\nTotal execution time: 49.2 seconds.\n\n# Predictions for missing (lost-to-follow-up) post-treatment values\n# Output format is a matrix with each column being a patient record and each row being one imputation\nimputations1 &lt;- predict(brmfit1, \n                        newdata=count_data %&gt;% \n                            filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n                            mutate(period=\"Post-treatment\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\n# Combine predictions with existing data frame of lost-to-follow-up records\n# (results in one row per imputation per patient), then bind together with \n# non-missing values (on-treatment and observed off-treatment values), for\n# which we have one record per patient and `.imputation` will be missing.\nimputed1 &lt;- count_data %&gt;% \n    filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations1[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n    dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period %in% c(\"On-treatment\", \"Post-treatment\") & \n                     follow_up&gt;0))\n\n\n\n7.5.2 Imputation based on the control group event rate (“jump-to-reference”)\nNow, we do an imputation based on control group data (i.e. assuming that after treatment discontinuation patients are off treatment and we use either the observed off-treatment data or impute missing data based on the control group event rate).\nFor this the code looks very similar to what we did in the previous subsection. We just slightly change our imputation model, and now predict based on setting the period covariate to \"On-treatment\" and the treatment covariate to \"0\" (=control group).\n\nbrmfit2 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + logBASE,\n               data=count_data %&gt;% \n                   filter(period==\"On-treatment\" & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 10.3 seconds.\nChain 2 finished in 10.2 seconds.\nChain 3 finished in 12.9 seconds.\nChain 4 finished in 12.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 11.4 seconds.\nTotal execution time: 46.3 seconds.\n\nimputations2 &lt;- predict(brmfit2, \n                        newdata=count_data %&gt;% \n                            filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n                        mutate(period=\"On-treatment\",\n                               treatment=\"0\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\nimputed2 &lt;- count_data %&gt;% \n    filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations2[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n  dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period %in% c(\"On-treatment\", \"Post-treatment\") & \n                     follow_up&gt;0))\n\n\n\n7.5.3 Imputation under a hypothetical estimand\nFinally, we also do an imputation based on on-treatment data of the same treatment group. I.e. we do an analysis under a hypothetical estimand “as if all patients had finished treatment”. As before, we just change our imputation model and the covariates for the data to be imputed accordingly.\n\nbrmfit3 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + logBASE,\n               data=count_data %&gt;% \n                   filter(period==\"On-treatment\" & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 12.6 seconds.\nChain 2 finished in 9.7 seconds.\nChain 3 finished in 10.4 seconds.\nChain 4 finished in 13.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 11.4 seconds.\nTotal execution time: 46.0 seconds.\n\nimputations3 &lt;- predict(brmfit3, \n                        newdata=count_data %&gt;% \n                            filter(period %in% c(\"Lost-to-follow-up\", \"Post-treatment\") & follow_up&gt;0) %&gt;%\n                         mutate(period=\"On-treatment\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\nimputed3 &lt;- count_data %&gt;% \n    filter(period %in% c(\"Lost-to-follow-up\", \"Post-treatment\") & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations3[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n    dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period==\"On-treatment\" &  follow_up&gt;0))\n\n\n\n7.5.4 Combining the imputations and a quick look at them\nWe combine the imputations from the different approaches and have a look at the mean rates under the different imputation approaches. As you can see, assuming a continued treatment effect results in a lower imputed mean event rate than assuming a cessation of the treatment effect (or basing the imputation on observed off-treatment data).\n\n# How different are the numbers?\nimputed &lt;- bind_rows(\n  imputed1 %&gt;%\n    mutate(Method=\"Retrieved drop-out imputation\"),\n  imputed2 %&gt;%\n    mutate(Method=\"Jump-to-reference (J2R)\"),\n  imputed3 %&gt;%\n    mutate(Method=\"Hypothetical (MAR)\")\n) \n\nimputed %&gt;%\n  group_by(Method, treatment, period, patient) %&gt;%\n  summarize(event_rate=mean(events/follow_up), .groups=\"drop\") %&gt;%\n  group_by(treatment, period, Method) %&gt;%\n  summarize(event_rate=mean(event_rate), .groups=\"drop\")\n\n# A tibble: 18 × 4\n   treatment period            Method                        event_rate\n   &lt;fct&gt;     &lt;fct&gt;             &lt;chr&gt;                              &lt;dbl&gt;\n 1 0         On-treatment      Hypothetical (MAR)                 1.61 \n 2 0         On-treatment      Jump-to-reference (J2R)            1.61 \n 3 0         On-treatment      Retrieved drop-out imputation      1.61 \n 4 0         Post-treatment    Hypothetical (MAR)                 1.68 \n 5 0         Post-treatment    Jump-to-reference (J2R)            1.48 \n 6 0         Post-treatment    Retrieved drop-out imputation      1.48 \n 7 0         Lost-to-follow-up Hypothetical (MAR)                 1.56 \n 8 0         Lost-to-follow-up Jump-to-reference (J2R)            1.55 \n 9 0         Lost-to-follow-up Retrieved drop-out imputation      1.54 \n10 1         On-treatment      Hypothetical (MAR)                 0.830\n11 1         On-treatment      Jump-to-reference (J2R)            0.830\n12 1         On-treatment      Retrieved drop-out imputation      0.830\n13 1         Post-treatment    Hypothetical (MAR)                 0.859\n14 1         Post-treatment    Jump-to-reference (J2R)            1.49 \n15 1         Post-treatment    Retrieved drop-out imputation      1.49 \n16 1         Lost-to-follow-up Hypothetical (MAR)                 0.899\n17 1         Lost-to-follow-up Jump-to-reference (J2R)            1.75 \n18 1         Lost-to-follow-up Retrieved drop-out imputation      1.49 \n\n\n\n\n7.5.5 Negative binomial regression for each imputed dataset\nWe can now conduct an analysis of each imputed dataset and look at how much the results vary between them. To do that, we need to fit a negative binomial regression to each imputed dataset. Negative binomial regression is implemented in all popular statistical software, but there are some particular challenges. Skip the next sub-section, if you are not interested in these particular points regarding negtive binomial regression using R.\n\n7.5.5.1 Negative binomial regression in R (optional subsection)\nFor example, the SAS/STATsoftware provides the GENMOD procedure and for the Poisson random effects version the COUNTREG procedure. Similarly, R has the glm.nb function in the MASS package. \nThese implementations differ in their parameterization. The GENMOD procedure estimates \\(\\kappa\\) and reverts to Poisson regression on the boundary of the parameter space when \\(\\hat{\\kappa}=0\\). The glm.nb function parameterizes the model in terms of \\(\\theta:=1/\\kappa\\) and its algorithm will not converge to a finite estimate \\(\\hat{\\theta}\\) whenever the maximum likelihood estimator for \\(\\kappa\\) is \\(\\hat{\\kappa}=0\\). However, after a sufficient number of iterations \\(\\hat{\\theta}\\) will be large enough that the estimates of the regression coefficients will be those we would obtain by Poisson regression — alternatively, we can use the results from Poisson regression in this case. This is implemented in the code below.\nInstead of the preferable observed Fisher information used in the GENMOD procedure, the glm.nb function uses the expected Fisher information for the standard errors of the regression coefficients — while the standard error for \\(\\hat{\\theta}\\) is based on the observed rather than the expected Fisher information — as of version 7.3-47 of the MASS package. Especially for small samples sizes this may result in too small standard errors. Thus, they should — especially for confirmatory clinical trials — be corrected, which is easily done using code provded by Bartlett.  \n\nglm.nb.cov &lt;- function(mod) {\n  # Basis of code for this function: \n  # https://stats.stackexchange.com/questions/221648/negative-binomial-regression-in-r-allowing-for-correlation-between-dispersion\n  #given a model fitted by glm.nb in MASS, this function returns a variance covariance matrix for the\n  #regression coefficients and dispersion parameter, without assuming independence between these\n  #note that the model must have been fitted with x=TRUE argument so that design matrix is available\n  #formulae based on p23-p24 of http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf\n  #and http://www.math.mcgill.ca/~dstephens/523/Papers/Lawless-1987-CJS.pdf\n  \n  k &lt;- mod$theta\n  p &lt;- dim(vcov(mod))[1]\n  \n  #construct observed information matrix\n  obsInfo &lt;- array(0, dim=c(p+1, p+1))\n  \n  #first calculate top left part for regression coefficients\n  for (i in 1:p) {\n    for (j in 1:p) {\n      obsInfo[i,j] &lt;- sum( (1+mod$y/mod$theta)*mod$fitted.values*mod$x[,i]*mod$x[,j] / (1+mod$fitted.values/mod$theta)^2  )\n    }\n  }\n  \n  #information for dispersion parameter\n  obsInfo[(p+1),(p+1)] &lt;- -sum(trigamma(mod$theta+mod$y) - trigamma(mod$theta) -\n                                 1/(mod$fitted.values+mod$theta) + (mod$theta+mod$y)/(mod$theta+mod$fitted.values)^2 - \n                                 1/(mod$fitted.values+mod$theta) + 1/mod$theta)\n  \n  #covariance between regression coefficients and dispersion\n  for (i in 1:p) {\n    obsInfo[(p+1),i] &lt;- -sum(((mod$y-mod$fitted.values) * mod$fitted.values / ( (mod$theta+mod$fitted.values)^2 )) * mod$x[,i] )\n    obsInfo[i,(p+1)] &lt;- obsInfo[(p+1),i]\n  }\n  \n  #return variance covariance matrix\n  solve(obsInfo)\n}\n\nfitmodel &lt;- function(data) {\n  # Function to fit NB model - returns result of Poisson model,\n  # if estimate of 1/dispersion parameter appears to go towards infinity.\n  # Additionally, standard errors are calculated using the observed\n  # information matrix by calling the previously defined function for that.\n  poifit &lt;- glm(\n    data = data,\n    formula = events ~ treatment + logBASE + offset(logfollowup),\n    control = glm.control(maxit = 2500),\n    family = poisson()\n  )\n  # glm.nb may fail to converge with a message such as\n  # \"Warning: step size truncated due to divergence\"\n  # \"Error in glm.fitter(X[, \"(Intercept)\", drop = FALSE], Y, w, offset = offset,: NA/NaN/Inf in 'x'\"\n  # We need to make sure we catch that issue. In that case we set nbfit = list(theta=10001), which\n  # results in Poisson regression results being used (see check on nbfit$theta later).\n  nbfit &lt;- tryCatch(\n      glm.nb(\n        data = data,\n        formula = aval ~ events ~ treatment + logBASE + offset(logfollowup),\n        start = poifit$coefficients,\n        init.theta = 1.25,\n        control = glm.control(maxit = 30),\n        x = TRUE\n      ), warning = function(err)\n        list(theta = 99999),\n      error = function(err)\n        list(theta = 99999)\n    )\n  \n  emm_options(msg.nesting = FALSE,\n              msg.interaction = FALSE)\n  if (nbfit$theta &gt; 10000) {\n    results &lt;- data.frame(\n      treatment = sort(unique(data$treatment)),\n      logrr = c(0, summary(poifit)$coefficients[2, 1]),\n      logrrSE = c(NA, summary(poifit)$coefficients[2, 2]),\n      as.matrix(data.frame(summary(\n        emmeans(\n          poifit,\n          ~ factor(treatment),\n          weights = \"proportional\",\n          at = list(logfollowup = 0)\n        )\n      ))[, c(\"emmean\", \"SE\")])\n    )\n  } else {\n    results &lt;- data.frame(\n      trt = sort(unique(data$treatment)),\n      logrr = c(0, summary(nbfit)$coefficients[2:2, 1]),\n      logrrSE = c(NA, sqrt(diag(\n        glm.nb.cov(nbfit)\n      ))[2:2]),\n      as.matrix(data.frame(summary(\n        emmeans(\n          nbfit,\n          ~ factor(treatment),\n          weights = \"proportional\",\n          at = list(logfollowup = 0),\n          vcov = glm.nb.cov(nbfit)[1:length(nbfit$coefficients), 1:length(nbfit$coefficients)]\n        )\n      ))[, c(\"emmean\", \"SE\")])\n    )\n    \n  }\n  \n  return(results)\n}\n\n\n\n7.5.5.2 Fitting the negative binomial regression to the imputed data\nActually fitting negative binomial regression on each imputed datasets is now easy with the function above:\n\nresults &lt;- tibble(.imputation = 1:max(imputed$.imputation, na.rm = T)) %&gt;%\n  mutate(results = parallel::mclapply(.imputation , function(x)\n    bind_rows(\n      imputed %&gt;%\n        filter(\n          Method == \"Hypothetical (MAR)\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Hypothetical (MAR)\"),\n      imputed %&gt;%\n        filter(\n          Method == \"Retrieved drop-out imputation\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Retrieved drop-out imputation\"),\n      imputed %&gt;%\n        filter(\n          Method == \"Jump-to-reference (J2R)\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Jump-to-reference (J2R)\")\n    ))) %&gt;%\n  unnest(results)\n\nresults %&gt;%\n  ggplot(aes(x=logrr, xmin=logrr-qnorm(0.975)*logrrSE,\n             xmax=logrr+qnorm(0.975)*logrrSE, y=.imputation, )) +\n  theme_bw(base_size=16) +\n  geom_vline(xintercept=0, linetype=2) +\n  geom_point() +\n  geom_errorbarh() +\n  facet_wrap(~Method) +\n  scale_x_continuous(limits = c(-0.9, 0)) +\n  ylab(\"Imputation number\") +\n  xlab(\"log-rate ratio (drug versus placebo, &lt;0 favors drug)\")\n\n\n\n\n\n\n\n\n\n\n\n7.5.6 Combining the results of the analysis of each dataset using Rubin’s rule\nFinally, we use Rubin’s rule to combine results across the different imputed datasets for each imputation method.\n\nrubins_rule &lt;- function(estimates, SEs) {\n  nmi &lt;- length(estimates)\n  mi_estimate &lt;- mean(estimates)\n  mi_SE &lt;- sqrt( sum(SEs^2)/nmi + (1+1/nmi) * sum((mi_estimate-estimates)^2) / (nmi-1) )\n  mi_Z &lt;- mi_estimate/mi_SE\n  mi_p &lt;- 2*min(pnorm(mi_Z), 1-pnorm(mi_Z))\n  return( tibble(logrr=mi_estimate, logrrSE=mi_SE, p=mi_p) )\n}\n\nfinal_results &lt;- results %&gt;%\n  group_by(Method) %&gt;%\n  do(rubins_rule(.$logrr, .$logrrSE)) %&gt;%\n  ungroup() %&gt;%\n  bind_rows(\n    count_data %&gt;%\n      filter(period==\"On-treatment\" & follow_up&gt;0) %&gt;%\n      mutate(logfollowup=log(follow_up)) %&gt;%\n      fitmodel(data=.) %&gt;% \n      as_tibble() %&gt;%\n      filter(treatment==\"1\") %&gt;% \n      mutate(Method=\"Hypothetical\\n(NegBin for on-treatment data)\",\n             p=2*min(pnorm(logrr/logrrSE), 1-pnorm(logrr/logrrSE)) ) %&gt;%\n      dplyr::select(-treatment, -emmean, -SE)\n  ) %&gt;%\n  mutate(upper = logrr + qnorm(0.975)*logrrSE,\n         lower = logrr - qnorm(0.975)*logrrSE)\n\n\n\n7.5.7 What if we want to do a Bayesian analysis after MI?\nIf we have multiple imputated datasets and wish to perform Bayesian inference, we would recommend the approach of described in the Bayesian Data Analysis book by Gelman et al. (Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B., 2014. Bayesian data analysis (3rd ed.). CRC press. p.452) that uses the (equally weighted) mixture distribution of the posterior distributions as the combined posterior distribution. I.e. we use all the (or a random subset) of the MCMC samples from each analysis as samples from the combined posterior distribution. This has been reported to perform well with a sufficient number of imputations such as 100 (Zhou, X. and Reiter, J.P., 2010. A note on Bayesian inference after multiple imputation. The American Statistician, 64(2), pp.159-163.).\nWe could of course write some code to fit separate Bayesian models for each imputed dataset and then combine the posterior samples, but we do not have to, because brms can take care of it for us. E.g. the code below would run 1 MCMC chain for each of 20 imputed datasets (without specifying any prior distributions, which we would in practice of course also do).\n\nlibrary(future)\nplan(multisession)\nfit_multiple1 &lt;- brm_multiple(\n    events | rate(follow_up) ~ treatment + logBASE,\n    data = map(1L:20L, \n               function(x) imputed1 %&gt;% \n                 filter(is.na(.imputation) | .imputation==x) %&gt;%\n                 group_by(patient, logBASE, treatment) %&gt;%\n                 summarize(events=sum(events), follow_up=sum(follow_up))),\n    family = negbinomial(),\n    chains = 1, seed = 6234)",
+    "text": "7.5 Implementation\n\n7.5.1 Imputation from observed off-treatment data\nFirst, let us conduct an imputation model based on the observed off-treatment data, in which we specify that event rates can be different on- and off-treatment for each treatment group. Using this model, we perform predictions for the missing post-treatment values (patients lost-to-follow-up).\nWe use family=poisson() here, because a Poisson model with a random patient effect on the intercept is approximately equivalent to a negative binomial regression model and this is a way of reflecting the correlation of multiple observation periods for a patient that is consistent with that. However, we could also make the imputation model more complex e.g. by allowing counts to be overdispersed within patient by using family=negbinomial(), or allowing for some interactions such as treatment:logBASE or period:logBASE.\n\n# Fit imputation model\nbrmfit1 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + period + treatment:period + logBASE,\n               data=count_data %&gt;% \n                   filter(period %in% c(\"On-treatment\", \"Post-treatment\") & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 11.2 seconds.\nChain 2 finished in 10.6 seconds.\nChain 3 finished in 13.5 seconds.\nChain 4 finished in 13.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 12.1 seconds.\nTotal execution time: 48.7 seconds.\n\n# Predictions for missing (lost-to-follow-up) post-treatment values\n# Output format is a matrix with each column being a patient record and each row being one imputation\nimputations1 &lt;- predict(brmfit1, \n                        newdata=count_data %&gt;% \n                            filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n                            mutate(period=\"Post-treatment\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\n# Combine predictions with existing data frame of lost-to-follow-up records\n# (results in one row per imputation per patient), then bind together with \n# non-missing values (on-treatment and observed off-treatment values), for\n# which we have one record per patient and `.imputation` will be missing.\nimputed1 &lt;- count_data %&gt;% \n    filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations1[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n    dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period %in% c(\"On-treatment\", \"Post-treatment\") & \n                     follow_up&gt;0))\n\n\n\n7.5.2 Imputation based on the control group event rate (“jump-to-reference”)\nNow, we do an imputation based on control group data (i.e. assuming that after treatment discontinuation patients are off treatment and we use either the observed off-treatment data or impute missing data based on the control group event rate).\nFor this the code looks very similar to what we did in the previous subsection. We just slightly change our imputation model, and now predict based on setting the period covariate to \"On-treatment\" and the treatment covariate to \"0\" (=control group).\n\nbrmfit2 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + logBASE,\n               data=count_data %&gt;% \n                   filter(period==\"On-treatment\" & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 10.1 seconds.\nChain 2 finished in 10.1 seconds.\nChain 3 finished in 12.8 seconds.\nChain 4 finished in 12.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 11.3 seconds.\nTotal execution time: 45.5 seconds.\n\nimputations2 &lt;- predict(brmfit2, \n                        newdata=count_data %&gt;% \n                            filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n                        mutate(period=\"On-treatment\",\n                               treatment=\"0\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\nimputed2 &lt;- count_data %&gt;% \n    filter(period %in% \"Lost-to-follow-up\" & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations2[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n  dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period %in% c(\"On-treatment\", \"Post-treatment\") & \n                     follow_up&gt;0))\n\n\n\n7.5.3 Imputation under a hypothetical estimand\nFinally, we also do an imputation based on on-treatment data of the same treatment group. I.e. we do an analysis under a hypothetical estimand “as if all patients had finished treatment”. As before, we just change our imputation model and the covariates for the data to be imputed accordingly.\n\nbrmfit3 &lt;- brm(formula = events | rate(follow_up) ~ 1 + (1|patient) + treatment + logBASE,\n               data=count_data %&gt;% \n                   filter(period==\"On-treatment\" & follow_up&gt;0),\n               family=poisson(),\n               refresh = 0, silent = 0)\n\nStart sampling\n\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 12.4 seconds.\nChain 2 finished in 9.6 seconds.\nChain 3 finished in 10.2 seconds.\nChain 4 finished in 12.8 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 11.3 seconds.\nTotal execution time: 45.4 seconds.\n\nimputations3 &lt;- predict(brmfit3, \n                        newdata=count_data %&gt;% \n                            filter(period %in% c(\"Lost-to-follow-up\", \"Post-treatment\") & follow_up&gt;0) %&gt;%\n                         mutate(period=\"On-treatment\"),\n                        summary = F,\n                        ndraws = num_imputations) \n\nimputed3 &lt;- count_data %&gt;% \n    filter(period %in% c(\"Lost-to-follow-up\", \"Post-treatment\") & follow_up&gt;0) %&gt;%\n    mutate(row_number=1:n(),\n           imputed_values = map(row_number, \n                                function(x) tibble(events=imputations3[,x]) %&gt;% \n                                            mutate(.imputation=1:n()))) %&gt;%\n    dplyr::select(-events) %&gt;%\n    unnest(imputed_values) %&gt;%\n    bind_rows(count_data %&gt;% \n              filter(period==\"On-treatment\" &  follow_up&gt;0))\n\n\n\n7.5.4 Combining the imputations and a quick look at them\nWe combine the imputations from the different approaches and have a look at the mean rates under the different imputation approaches. As you can see, assuming a continued treatment effect results in a lower imputed mean event rate than assuming a cessation of the treatment effect (or basing the imputation on observed off-treatment data).\n\n# How different are the numbers?\nimputed &lt;- bind_rows(\n  imputed1 %&gt;%\n    mutate(Method=\"Retrieved drop-out imputation\"),\n  imputed2 %&gt;%\n    mutate(Method=\"Jump-to-reference (J2R)\"),\n  imputed3 %&gt;%\n    mutate(Method=\"Hypothetical (MAR)\")\n) \n\nimputed %&gt;%\n  group_by(Method, treatment, period, patient) %&gt;%\n  summarize(event_rate=mean(events/follow_up), .groups=\"drop\") %&gt;%\n  group_by(treatment, period, Method) %&gt;%\n  summarize(event_rate=mean(event_rate), .groups=\"drop\")\n\n# A tibble: 18 × 4\n   treatment period            Method                        event_rate\n   &lt;fct&gt;     &lt;fct&gt;             &lt;chr&gt;                              &lt;dbl&gt;\n 1 0         On-treatment      Hypothetical (MAR)                 1.61 \n 2 0         On-treatment      Jump-to-reference (J2R)            1.61 \n 3 0         On-treatment      Retrieved drop-out imputation      1.61 \n 4 0         Post-treatment    Hypothetical (MAR)                 1.68 \n 5 0         Post-treatment    Jump-to-reference (J2R)            1.48 \n 6 0         Post-treatment    Retrieved drop-out imputation      1.48 \n 7 0         Lost-to-follow-up Hypothetical (MAR)                 1.56 \n 8 0         Lost-to-follow-up Jump-to-reference (J2R)            1.55 \n 9 0         Lost-to-follow-up Retrieved drop-out imputation      1.54 \n10 1         On-treatment      Hypothetical (MAR)                 0.830\n11 1         On-treatment      Jump-to-reference (J2R)            0.830\n12 1         On-treatment      Retrieved drop-out imputation      0.830\n13 1         Post-treatment    Hypothetical (MAR)                 0.859\n14 1         Post-treatment    Jump-to-reference (J2R)            1.49 \n15 1         Post-treatment    Retrieved drop-out imputation      1.49 \n16 1         Lost-to-follow-up Hypothetical (MAR)                 0.899\n17 1         Lost-to-follow-up Jump-to-reference (J2R)            1.75 \n18 1         Lost-to-follow-up Retrieved drop-out imputation      1.49 \n\n\n\n\n7.5.5 Negative binomial regression for each imputed dataset\nWe can now conduct an analysis of each imputed dataset and look at how much the results vary between them. To do that, we need to fit a negative binomial regression to each imputed dataset. Negative binomial regression is implemented in all popular statistical software, but there are some particular challenges. Skip the next sub-section, if you are not interested in these particular points regarding negtive binomial regression using R.\n\n7.5.5.1 Negative binomial regression in R (optional subsection)\nFor example, the SAS/STATsoftware provides the GENMOD procedure and for the Poisson random effects version the COUNTREG procedure. Similarly, R has the glm.nb function in the MASS package. \nThese implementations differ in their parameterization. The GENMOD procedure estimates \\(\\kappa\\) and reverts to Poisson regression on the boundary of the parameter space when \\(\\hat{\\kappa}=0\\). The glm.nb function parameterizes the model in terms of \\(\\theta:=1/\\kappa\\) and its algorithm will not converge to a finite estimate \\(\\hat{\\theta}\\) whenever the maximum likelihood estimator for \\(\\kappa\\) is \\(\\hat{\\kappa}=0\\). However, after a sufficient number of iterations \\(\\hat{\\theta}\\) will be large enough that the estimates of the regression coefficients will be those we would obtain by Poisson regression — alternatively, we can use the results from Poisson regression in this case. This is implemented in the code below.\nInstead of the preferable observed Fisher information used in the GENMOD procedure, the glm.nb function uses the expected Fisher information for the standard errors of the regression coefficients — while the standard error for \\(\\hat{\\theta}\\) is based on the observed rather than the expected Fisher information — as of version 7.3-47 of the MASS package. Especially for small samples sizes this may result in too small standard errors. Thus, they should — especially for confirmatory clinical trials — be corrected, which is easily done using code provded by Bartlett.  \n\nglm.nb.cov &lt;- function(mod) {\n  # Basis of code for this function: \n  # https://stats.stackexchange.com/questions/221648/negative-binomial-regression-in-r-allowing-for-correlation-between-dispersion\n  #given a model fitted by glm.nb in MASS, this function returns a variance covariance matrix for the\n  #regression coefficients and dispersion parameter, without assuming independence between these\n  #note that the model must have been fitted with x=TRUE argument so that design matrix is available\n  #formulae based on p23-p24 of http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf\n  #and http://www.math.mcgill.ca/~dstephens/523/Papers/Lawless-1987-CJS.pdf\n  \n  k &lt;- mod$theta\n  p &lt;- dim(vcov(mod))[1]\n  \n  #construct observed information matrix\n  obsInfo &lt;- array(0, dim=c(p+1, p+1))\n  \n  #first calculate top left part for regression coefficients\n  for (i in 1:p) {\n    for (j in 1:p) {\n      obsInfo[i,j] &lt;- sum( (1+mod$y/mod$theta)*mod$fitted.values*mod$x[,i]*mod$x[,j] / (1+mod$fitted.values/mod$theta)^2  )\n    }\n  }\n  \n  #information for dispersion parameter\n  obsInfo[(p+1),(p+1)] &lt;- -sum(trigamma(mod$theta+mod$y) - trigamma(mod$theta) -\n                                 1/(mod$fitted.values+mod$theta) + (mod$theta+mod$y)/(mod$theta+mod$fitted.values)^2 - \n                                 1/(mod$fitted.values+mod$theta) + 1/mod$theta)\n  \n  #covariance between regression coefficients and dispersion\n  for (i in 1:p) {\n    obsInfo[(p+1),i] &lt;- -sum(((mod$y-mod$fitted.values) * mod$fitted.values / ( (mod$theta+mod$fitted.values)^2 )) * mod$x[,i] )\n    obsInfo[i,(p+1)] &lt;- obsInfo[(p+1),i]\n  }\n  \n  #return variance covariance matrix\n  solve(obsInfo)\n}\n\nfitmodel &lt;- function(data) {\n  # Function to fit NB model - returns result of Poisson model,\n  # if estimate of 1/dispersion parameter appears to go towards infinity.\n  # Additionally, standard errors are calculated using the observed\n  # information matrix by calling the previously defined function for that.\n  poifit &lt;- glm(\n    data = data,\n    formula = events ~ treatment + logBASE + offset(logfollowup),\n    control = glm.control(maxit = 2500),\n    family = poisson()\n  )\n  # glm.nb may fail to converge with a message such as\n  # \"Warning: step size truncated due to divergence\"\n  # \"Error in glm.fitter(X[, \"(Intercept)\", drop = FALSE], Y, w, offset = offset,: NA/NaN/Inf in 'x'\"\n  # We need to make sure we catch that issue. In that case we set nbfit = list(theta=10001), which\n  # results in Poisson regression results being used (see check on nbfit$theta later).\n  nbfit &lt;- tryCatch(\n      glm.nb(\n        data = data,\n        formula = aval ~ events ~ treatment + logBASE + offset(logfollowup),\n        start = poifit$coefficients,\n        init.theta = 1.25,\n        control = glm.control(maxit = 30),\n        x = TRUE\n      ), warning = function(err)\n        list(theta = 99999),\n      error = function(err)\n        list(theta = 99999)\n    )\n  \n  emm_options(msg.nesting = FALSE,\n              msg.interaction = FALSE)\n  if (nbfit$theta &gt; 10000) {\n    results &lt;- data.frame(\n      treatment = sort(unique(data$treatment)),\n      logrr = c(0, summary(poifit)$coefficients[2, 1]),\n      logrrSE = c(NA, summary(poifit)$coefficients[2, 2]),\n      as.matrix(data.frame(summary(\n        emmeans(\n          poifit,\n          ~ factor(treatment),\n          weights = \"proportional\",\n          at = list(logfollowup = 0)\n        )\n      ))[, c(\"emmean\", \"SE\")])\n    )\n  } else {\n    results &lt;- data.frame(\n      trt = sort(unique(data$treatment)),\n      logrr = c(0, summary(nbfit)$coefficients[2:2, 1]),\n      logrrSE = c(NA, sqrt(diag(\n        glm.nb.cov(nbfit)\n      ))[2:2]),\n      as.matrix(data.frame(summary(\n        emmeans(\n          nbfit,\n          ~ factor(treatment),\n          weights = \"proportional\",\n          at = list(logfollowup = 0),\n          vcov = glm.nb.cov(nbfit)[1:length(nbfit$coefficients), 1:length(nbfit$coefficients)]\n        )\n      ))[, c(\"emmean\", \"SE\")])\n    )\n    \n  }\n  \n  return(results)\n}\n\n\n\n7.5.5.2 Fitting the negative binomial regression to the imputed data\nActually fitting negative binomial regression on each imputed datasets is now easy with the function above:\n\nresults &lt;- tibble(.imputation = 1:max(imputed$.imputation, na.rm = T)) %&gt;%\n  mutate(results = parallel::mclapply(.imputation , function(x)\n    bind_rows(\n      imputed %&gt;%\n        filter(\n          Method == \"Hypothetical (MAR)\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Hypothetical (MAR)\"),\n      imputed %&gt;%\n        filter(\n          Method == \"Retrieved drop-out imputation\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Retrieved drop-out imputation\"),\n      imputed %&gt;%\n        filter(\n          Method == \"Jump-to-reference (J2R)\" &\n            .imputation %in% c(x, NA_integer_)\n        ) %&gt;%\n        group_by(patient, treatment, BASE, logBASE, .imputation) %&gt;%\n        summarize(\n          logfollowup = log(sum(follow_up)),\n          events = sum(events),\n          .groups = \"drop\"\n        ) %&gt;%\n        fitmodel(data = .) %&gt;%\n        as_tibble() %&gt;%\n        filter(treatment == \"1\") %&gt;%\n        mutate(Method = \"Jump-to-reference (J2R)\")\n    ))) %&gt;%\n  unnest(results)\n\nresults %&gt;%\n  ggplot(aes(x=logrr, xmin=logrr-qnorm(0.975)*logrrSE,\n             xmax=logrr+qnorm(0.975)*logrrSE, y=.imputation, )) +\n  theme_bw(base_size=16) +\n  geom_vline(xintercept=0, linetype=2) +\n  geom_point() +\n  geom_errorbarh() +\n  facet_wrap(~Method) +\n  scale_x_continuous(limits = c(-0.9, 0)) +\n  ylab(\"Imputation number\") +\n  xlab(\"log-rate ratio (drug versus placebo, &lt;0 favors drug)\")\n\n\n\n\n\n\n\n\n\n\n\n7.5.6 Combining the results of the analysis of each dataset using Rubin’s rule\nFinally, we use Rubin’s rule to combine results across the different imputed datasets for each imputation method.\n\nrubins_rule &lt;- function(estimates, SEs) {\n  nmi &lt;- length(estimates)\n  mi_estimate &lt;- mean(estimates)\n  mi_SE &lt;- sqrt( sum(SEs^2)/nmi + (1+1/nmi) * sum((mi_estimate-estimates)^2) / (nmi-1) )\n  mi_Z &lt;- mi_estimate/mi_SE\n  mi_p &lt;- 2*min(pnorm(mi_Z), 1-pnorm(mi_Z))\n  return( tibble(logrr=mi_estimate, logrrSE=mi_SE, p=mi_p) )\n}\n\nfinal_results &lt;- results %&gt;%\n  group_by(Method) %&gt;%\n  do(rubins_rule(.$logrr, .$logrrSE)) %&gt;%\n  ungroup() %&gt;%\n  bind_rows(\n    count_data %&gt;%\n      filter(period==\"On-treatment\" & follow_up&gt;0) %&gt;%\n      mutate(logfollowup=log(follow_up)) %&gt;%\n      fitmodel(data=.) %&gt;% \n      as_tibble() %&gt;%\n      filter(treatment==\"1\") %&gt;% \n      mutate(Method=\"Hypothetical\\n(NegBin for on-treatment data)\",\n             p=2*min(pnorm(logrr/logrrSE), 1-pnorm(logrr/logrrSE)) ) %&gt;%\n      dplyr::select(-treatment, -emmean, -SE)\n  ) %&gt;%\n  mutate(upper = logrr + qnorm(0.975)*logrrSE,\n         lower = logrr - qnorm(0.975)*logrrSE)\n\n\n\n7.5.7 What if we want to do a Bayesian analysis after MI?\nIf we have multiple imputated datasets and wish to perform Bayesian inference, we would recommend the approach of described in the Bayesian Data Analysis book by Gelman et al. (Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B., 2014. Bayesian data analysis (3rd ed.). CRC press. p.452) that uses the (equally weighted) mixture distribution of the posterior distributions as the combined posterior distribution. I.e. we use all the (or a random subset) of the MCMC samples from each analysis as samples from the combined posterior distribution. This has been reported to perform well with a sufficient number of imputations such as 100 (Zhou, X. and Reiter, J.P., 2010. A note on Bayesian inference after multiple imputation. The American Statistician, 64(2), pp.159-163.).\nWe could of course write some code to fit separate Bayesian models for each imputed dataset and then combine the posterior samples, but we do not have to, because brms can take care of it for us. E.g. the code below would run 1 MCMC chain for each of 20 imputed datasets (without specifying any prior distributions, which we would in practice of course also do).\n\nlibrary(future)\nplan(multisession)\nfit_multiple1 &lt;- brm_multiple(\n    events | rate(follow_up) ~ treatment + logBASE,\n    data = map(1L:20L, \n               function(x) imputed1 %&gt;% \n                 filter(is.na(.imputation) | .imputation==x) %&gt;%\n                 group_by(patient, logBASE, treatment) %&gt;%\n                 summarize(events=sum(events), follow_up=sum(follow_up))),\n    family = negbinomial(),\n    chains = 1, seed = 6234)",
     "crumbs": [
       "Case studies",
       "<span class='chapter-number'>7</span>  <span class='chapter-title'>Multiple imputation</span>"
@@ -703,7 +703,7 @@
     "href": "src/02i_time_to_event.html#data",
     "title": "9  Time-to-event data",
     "section": "9.4 Data",
-    "text": "9.4 Data\nWe demonstrate here the trial analysis using a fake data simulation of the trial design. As key modeling choices for the parametric modeling approach are motivated from the historical data we start with the presentation of the historical data set and then describe the details of the trial simulation.\n\n9.4.1 Historical data\nThe CheckMate 649 trial (Janjigian et al. 2021) included 789 gastric cancer patients which were treated with either chemoA or chemoB and their progression free survival (PFS) was reported. Both chemotherapies were used in a roughly 1:1 ratio such that we will consider the reported data as “average” between both chemo therapies (despite lack of randomization occured wrt to chemotherapy assignment).\nThe data of trial has here been reconstructed from the published survival curves:\n\n\n\n\n\n\n\n\n\nTo now evaluate if the Weibull probability density function is an appropiate choice for this data we consider if the non-parametric estimate of the survivor function is compatible with properties of Weibull distributed event times. Specifically, if we transform the survivor function \\(S_{\\mbox{Weibull}}(t)\\) of the Weibull probability density function with the complementary log-log transformation, we obtain that\n\\[ \\mbox{cloglog}(S_{\\mbox{Weibull}}(t)) = \\log( - \\log(S_{\\mbox{Weibull}}(t))) =\n-\\alpha \\, \\log(s)  + \\alpha \\, \\log(t)\\]\nholds. Therefore, we can visualize an estimate of the survivor function on a transformed scale as a function of \\(\\log(t)\\) and we should observe a straight line with slope \\(\\alpha\\) and intercept \\(-\\alpha \\, \\log(s)\\):\n\nkm &lt;- survfit(Surv(time, status) ~ 1, data=hdata2)\nggsurvplot(km, data=hdata2, fun = \"cloglog\")\n\n\n\n\n\n\n\n\nWhile the line is not perfectly straight over the entire follow-up, the assumption of a straight line within the time period including the bulk of all events falls in the range of a good approximation with a line (as can be seen from the numbers at risk of the plot above).\nAs an additional evaluation of the Weibull distributional form we model the CheckMate 649 data via brms and perform a posterior predictive check. The priors are specified as illustrated in the later section on the used priors. The result looks reasonable.\n\nmodel_weibull_hist &lt;- bf(time | cens(1-status) ~ 1, family=weibull())\n\nprior_weibull_hist &lt;- prior(normal(log(6), log(4)/1.64), class=Intercept) +\n    prior(gamma(3, 2.7), class=shape)\n\nfit_hist_checkmate &lt;- brm(model_weibull_hist,\n                          data=hdata2,\n                          prior=prior_weibull_hist,\n                          seed=234235,\n                          refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.5 seconds.\nChain 2 finished in 1.4 seconds.\nChain 3 finished in 1.4 seconds.\nChain 4 finished in 1.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.4 seconds.\nTotal execution time: 5.9 seconds.\n\nfit_hist_checkmate\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: time | cens(1 - status) ~ 1 \n   Data: hdata2 (Number of observations: 789) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept     2.44      0.04     2.37     2.51 1.00     2794     2455\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.22      0.04     1.14     1.30 1.00     2865     2574\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nWith the posterior we can now perform a posterior predictive check. We do so by drawing 100 samples of the model parameters from the posterior and draw for each sample a realization of the event times for all patients. Hence, for each sample one obtains a fake data set simulated according to the Weibull distribution. Finally, each fake data set is summarized by a Kaplan-Meier estimate and drawn in light color while the actually observed data is summarized in the same way and drawn using a thick line. If the fitted model is appropiate, then the observed Kaplan Meier estimate should look just like a random realization among the many fake data generated ones. However, we observe that over the entire follow-up time the observed data is leaving for late follow-up times the bulk of the fake data generated Kaplan Meier curves. This may suggest a deviation from the Weibull distribution or the possibility that time-varying covariates are helpful. Nonetheless, if one focusses on the follow-up time frame, which includes the majoritiy of events (up to 15 month), then the posterior predictive check does look reasonable. This is why we use the Weibull distribution without a time-varying covariate for simplicity.\n\np_full_fup &lt;- pp_check(fit_hist_checkmate, status_y=hdata2$status,\n             type=\"km_overlay\", ndraws=100) +\n    scale_y_continuous(breaks=seq(0,1,by=0.1)) +\n  xlab(\"Time [month]\") + coord_cartesian(xlim=c(0, 35))\np_sub_fup &lt;- p_full_fup + coord_cartesian(xlim=c(0, 15))\nggpubr::ggarrange(p_full_fup, p_sub_fup, common.legend = TRUE)\n\n\n\n\n\n\n\n\nNote that the above predictive check does not simulate the censoring process. Thus, the censoring process is ignored and non-informative censoring is thereby assumed.\n\n\n9.4.2 Trial simulation\nHere we consider a fake data simulation of the trial of interest evaluating a novel treatment of gastric cancer for which an established treatment is available already. The trial is a randomized trial to study the efficacy and safety of adding a drug CompoundX in combination with a monoclonal antibody mAb + standard of care (SoC). Two chemotherapy backbones, ChemoA and ChemoB were used as SoC. There are in total four treatment arms including two control groups, which are mAb+ChemoA (controlChemoA) and mAb+ChemoB (controlChemoB), and two corresponding active arms with the test treatment administrated in addition, that is CompoundX+mAb+ChemoA (activeChemoA), CompoundX+mAb+ChemoB (activeChemoB).\nSynthetic data were generated for the randomised study. We let the sample size of both active and control group be 100, respectively. The sample size has been chosen to ensure adequate trial power of 80%. Efficacy was assessed using progression-free survival (PFS) measured in units of months. For patients not experiencing a progression event at the time of the analysis, their event time will be right censored at the time of their last valid tumor assessment. For patients experiencing an event, their event time will be interval censored, with the upper limit of the interval being the study day on which the progression event is observed, and the lower limit being the day of the last preceding valid tumor assessment at which the patient’s progression-free status was confirmed. While interval censoring is a supported feature of brms, this case study will for simplicity use the actual event times and leave it to the reader to explore interval censoring with brms.\nThere are two historical datasets available. The first one with 400 patients on a control arm, mAb+ChemoA (controlChemoA), is simulated together with the randomised study. And the PFS data of 789 patients receiving first-line programmed cell death (PD-1) inhibitor+chemotherapy are available from a randomised, open-label, phase 3 trial, CheckMate 649 (Janjigian et al. 2021), and it will be introduced later. These data are highly relevant and we would like to integrate the historical information into our analysis.\nSimulated data will be generated using simsurv function via simulating survival times from the Weibull model using the proportional hazard form of the model. We assumed a 20% hazard rate (HR) reduction due to active treatment administration (CompoundX) to the active group in comparison to the control group and a 5% difference in HR for the two chemotherapy SoCs (ChemoB is 5% worse than ChemoA on the hazard scale). We also assumed that both of the historical data have a 2% worse outcome on the hazards scale and we represent the two studies with hist1 and hist2 respectively (hist1 for simulated historidal data and hist2 for CheckMate 649).\n\n\nShow the code\nset.seed(46767)\n# n per group\nn_grp &lt;- 100\nn_hist &lt;- 400\n# use month as time-unit\nrate_1 &lt;- 1 / 6\nrate_cens &lt;- 1 / 10\nbeta &lt;- c(trt=-0.2,        ## roughly 20% HR reduction\n          soc_alt=0.05,    ## alternative chemotherapy is 5% worse\n          hist1=0.02)      ## simulated historical data has a 2% worse outcome\n# covariates of simulated trial data\ncovs &lt;- data.frame(id = seq(1, 2*n_grp), trt = c(0L, 1L),\n                   soc_alt=rbinom(2*n_grp, 1L, 0.3), hist1=0L, hist2=0L)\n\n# covariates of historical data\nhcovs &lt;- data.frame(id = seq(2*n_grp+1, 2*n_grp+1 +n_hist - 1),\n                    trt = c(0), soc_alt=0, hist1=1L, hist2=0L)\nsimulate_trial &lt;- function(lambda, gamma, lambda_cens, x, betas) {\n    ## simulate censoring times, note that we do not simulate end of\n    ## trial with maxt for now. Also note the different parametrization of\n    ## simsurv which models log(shape) with a linear predictor.\n    cens &lt;- simsurv(lambdas = lambda_cens, gammas = 1, x=x)\n    events &lt;- simsurv(lambdas = lambda, gammas = gamma, x=x, betas=betas)\n    names(cens) &lt;- paste0(names(cens), \"_cens\")\n    bind_cols(events, select(cens, -id_cens), select(x, -id)) %&gt;%\n        rename(censtime=eventtime_cens) %&gt;%\n    mutate(event=1L*(eventtime &lt;= censtime),\n           y=if_else(event==1, eventtime, censtime),\n           status=NULL, status_cens=NULL) %&gt;%\n        relocate(id, y, event) %&gt;%\n        mutate(soc=factor(soc_alt, c(0, 1), c(\"ChemoA\", \"ChemoB\")),\n               trt_ind=trt,\n               trt=factor(trt_ind, c(0,1), c(\"control\", \"active\")),\n               arm=factor(paste0(c(\"control\", \"active\")[trt_ind + 1], soc),\n                          levels=c(\"activeChemoA\", \"controlChemoA\",\n                                   \"activeChemoB\", \"controlChemoB\")))\n}\n\n# using shape=1 for simulation (corresponding to exponential survival times)\nsim &lt;- simulate_trial(rate_1, 1, rate_cens, covs, beta)\nhdata1 &lt;- simulate_trial(rate_1, 1, rate_cens, hcovs, beta) %&gt;%\n  mutate(id=id+max(sim$id))\n\n\nWe can check how the simulated data looks like and its sample size.\n\ngt_preview(sim) %&gt;% fmt_number(where(is.double), decimals=1)\n\n\n\n\n\n\n\n\n\nid\ny\nevent\neventtime\ncenstime\ntrt\nsoc_alt\nhist1\nhist2\nsoc\ntrt_ind\narm\n\n\n\n\n1\n1\n7.7\n0\n12.3\n7.7\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n2\n2\n0.1\n0\n20.4\n0.1\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n3\n3\n4.7\n0\n10.3\n4.7\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n4\n4\n2.7\n0\n19.0\n2.7\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n5\n5\n3.6\n1\n3.6\n7.5\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n6..199\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n200\n200\n11.8\n0\n17.1\n11.8\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n\n\n\n\n\ntable(sim$arm)\n\n\n activeChemoA controlChemoA  activeChemoB controlChemoB \n           72            73            28            27 \n\ntable(hdata1$arm)\n\n\n activeChemoA controlChemoA  activeChemoB controlChemoB \n            0           400             0             0",
+    "text": "9.4 Data\nWe demonstrate here the trial analysis using a fake data simulation of the trial design. As key modeling choices for the parametric modeling approach are motivated from the historical data we start with the presentation of the historical data set and then describe the details of the trial simulation.\n\n9.4.1 Historical data\nThe CheckMate 649 trial (Janjigian et al. 2021) included 789 gastric cancer patients which were treated with either chemoA or chemoB and their progression free survival (PFS) was reported. Both chemotherapies were used in a roughly 1:1 ratio such that we will consider the reported data as “average” between both chemo therapies (despite lack of randomization occured wrt to chemotherapy assignment).\nThe data of trial has here been reconstructed from the published survival curves:\n\n\n\n\n\n\n\n\n\nTo now evaluate if the Weibull probability density function is an appropiate choice for this data we consider if the non-parametric estimate of the survivor function is compatible with properties of Weibull distributed event times. Specifically, if we transform the survivor function \\(S_{\\mbox{Weibull}}(t)\\) of the Weibull probability density function with the complementary log-log transformation, we obtain that\n\\[ \\mbox{cloglog}(S_{\\mbox{Weibull}}(t)) = \\log( - \\log(S_{\\mbox{Weibull}}(t))) =\n-\\alpha \\, \\log(s)  + \\alpha \\, \\log(t)\\]\nholds. Therefore, we can visualize an estimate of the survivor function on a transformed scale as a function of \\(\\log(t)\\) and we should observe a straight line with slope \\(\\alpha\\) and intercept \\(-\\alpha \\, \\log(s)\\):\n\nkm &lt;- survfit(Surv(time, status) ~ 1, data=hdata2)\nggsurvplot(km, data=hdata2, fun = \"cloglog\")\n\n\n\n\n\n\n\n\nWhile the line is not perfectly straight over the entire follow-up, the assumption of a straight line within the time period including the bulk of all events falls in the range of a good approximation with a line (as can be seen from the numbers at risk of the plot above).\nAs an additional evaluation of the Weibull distributional form we model the CheckMate 649 data via brms and perform a posterior predictive check. The priors are specified as illustrated in the later section on the used priors. The result looks reasonable.\n\nmodel_weibull_hist &lt;- bf(time | cens(1-status) ~ 1, family=weibull())\n\nprior_weibull_hist &lt;- prior(normal(log(6), log(4)/1.64), class=Intercept) +\n    prior(gamma(3, 2.7), class=shape)\n\nfit_hist_checkmate &lt;- brm(model_weibull_hist,\n                          data=hdata2,\n                          prior=prior_weibull_hist,\n                          seed=234235,\n                          refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.3 seconds.\nChain 2 finished in 1.4 seconds.\nChain 3 finished in 1.4 seconds.\nChain 4 finished in 1.3 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.3 seconds.\nTotal execution time: 5.7 seconds.\n\nfit_hist_checkmate\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: time | cens(1 - status) ~ 1 \n   Data: hdata2 (Number of observations: 789) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept     2.44      0.04     2.37     2.51 1.00     2794     2455\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.22      0.04     1.14     1.30 1.00     2865     2574\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nWith the posterior we can now perform a posterior predictive check. We do so by drawing 100 samples of the model parameters from the posterior and draw for each sample a realization of the event times for all patients. Hence, for each sample one obtains a fake data set simulated according to the Weibull distribution. Finally, each fake data set is summarized by a Kaplan-Meier estimate and drawn in light color while the actually observed data is summarized in the same way and drawn using a thick line. If the fitted model is appropiate, then the observed Kaplan Meier estimate should look just like a random realization among the many fake data generated ones. However, we observe that over the entire follow-up time the observed data is leaving for late follow-up times the bulk of the fake data generated Kaplan Meier curves. This may suggest a deviation from the Weibull distribution or the possibility that time-varying covariates are helpful. Nonetheless, if one focusses on the follow-up time frame, which includes the majoritiy of events (up to 15 month), then the posterior predictive check does look reasonable. This is why we use the Weibull distribution without a time-varying covariate for simplicity.\n\np_full_fup &lt;- pp_check(fit_hist_checkmate, status_y=hdata2$status,\n             type=\"km_overlay\", ndraws=100) +\n    scale_y_continuous(breaks=seq(0,1,by=0.1)) +\n  xlab(\"Time [month]\") + coord_cartesian(xlim=c(0, 35))\np_sub_fup &lt;- p_full_fup + coord_cartesian(xlim=c(0, 15))\nggpubr::ggarrange(p_full_fup, p_sub_fup, common.legend = TRUE)\n\n\n\n\n\n\n\n\nNote that the above predictive check does not simulate the censoring process. Thus, the censoring process is ignored and non-informative censoring is thereby assumed.\n\n\n9.4.2 Trial simulation\nHere we consider a fake data simulation of the trial of interest evaluating a novel treatment of gastric cancer for which an established treatment is available already. The trial is a randomized trial to study the efficacy and safety of adding a drug CompoundX in combination with a monoclonal antibody mAb + standard of care (SoC). Two chemotherapy backbones, ChemoA and ChemoB were used as SoC. There are in total four treatment arms including two control groups, which are mAb+ChemoA (controlChemoA) and mAb+ChemoB (controlChemoB), and two corresponding active arms with the test treatment administrated in addition, that is CompoundX+mAb+ChemoA (activeChemoA), CompoundX+mAb+ChemoB (activeChemoB).\nSynthetic data were generated for the randomised study. We let the sample size of both active and control group be 100, respectively. The sample size has been chosen to ensure adequate trial power of 80%. Efficacy was assessed using progression-free survival (PFS) measured in units of months. For patients not experiencing a progression event at the time of the analysis, their event time will be right censored at the time of their last valid tumor assessment. For patients experiencing an event, their event time will be interval censored, with the upper limit of the interval being the study day on which the progression event is observed, and the lower limit being the day of the last preceding valid tumor assessment at which the patient’s progression-free status was confirmed. While interval censoring is a supported feature of brms, this case study will for simplicity use the actual event times and leave it to the reader to explore interval censoring with brms.\nThere are two historical datasets available. The first one with 400 patients on a control arm, mAb+ChemoA (controlChemoA), is simulated together with the randomised study. And the PFS data of 789 patients receiving first-line programmed cell death (PD-1) inhibitor+chemotherapy are available from a randomised, open-label, phase 3 trial, CheckMate 649 (Janjigian et al. 2021), and it will be introduced later. These data are highly relevant and we would like to integrate the historical information into our analysis.\nSimulated data will be generated using simsurv function via simulating survival times from the Weibull model using the proportional hazard form of the model. We assumed a 20% hazard rate (HR) reduction due to active treatment administration (CompoundX) to the active group in comparison to the control group and a 5% difference in HR for the two chemotherapy SoCs (ChemoB is 5% worse than ChemoA on the hazard scale). We also assumed that both of the historical data have a 2% worse outcome on the hazards scale and we represent the two studies with hist1 and hist2 respectively (hist1 for simulated historidal data and hist2 for CheckMate 649).\n\n\nShow the code\nset.seed(46767)\n# n per group\nn_grp &lt;- 100\nn_hist &lt;- 400\n# use month as time-unit\nrate_1 &lt;- 1 / 6\nrate_cens &lt;- 1 / 10\nbeta &lt;- c(trt=-0.2,        ## roughly 20% HR reduction\n          soc_alt=0.05,    ## alternative chemotherapy is 5% worse\n          hist1=0.02)      ## simulated historical data has a 2% worse outcome\n# covariates of simulated trial data\ncovs &lt;- data.frame(id = seq(1, 2*n_grp), trt = c(0L, 1L),\n                   soc_alt=rbinom(2*n_grp, 1L, 0.3), hist1=0L, hist2=0L)\n\n# covariates of historical data\nhcovs &lt;- data.frame(id = seq(2*n_grp+1, 2*n_grp+1 +n_hist - 1),\n                    trt = c(0), soc_alt=0, hist1=1L, hist2=0L)\nsimulate_trial &lt;- function(lambda, gamma, lambda_cens, x, betas) {\n    ## simulate censoring times, note that we do not simulate end of\n    ## trial with maxt for now. Also note the different parametrization of\n    ## simsurv which models log(shape) with a linear predictor.\n    cens &lt;- simsurv(lambdas = lambda_cens, gammas = 1, x=x)\n    events &lt;- simsurv(lambdas = lambda, gammas = gamma, x=x, betas=betas)\n    names(cens) &lt;- paste0(names(cens), \"_cens\")\n    bind_cols(events, select(cens, -id_cens), select(x, -id)) %&gt;%\n        rename(censtime=eventtime_cens) %&gt;%\n    mutate(event=1L*(eventtime &lt;= censtime),\n           y=if_else(event==1, eventtime, censtime),\n           status=NULL, status_cens=NULL) %&gt;%\n        relocate(id, y, event) %&gt;%\n        mutate(soc=factor(soc_alt, c(0, 1), c(\"ChemoA\", \"ChemoB\")),\n               trt_ind=trt,\n               trt=factor(trt_ind, c(0,1), c(\"control\", \"active\")),\n               arm=factor(paste0(c(\"control\", \"active\")[trt_ind + 1], soc),\n                          levels=c(\"activeChemoA\", \"controlChemoA\",\n                                   \"activeChemoB\", \"controlChemoB\")))\n}\n\n# using shape=1 for simulation (corresponding to exponential survival times)\nsim &lt;- simulate_trial(rate_1, 1, rate_cens, covs, beta)\nhdata1 &lt;- simulate_trial(rate_1, 1, rate_cens, hcovs, beta) %&gt;%\n  mutate(id=id+max(sim$id))\n\n\nWe can check how the simulated data looks like and its sample size.\n\ngt_preview(sim) %&gt;% fmt_number(where(is.double), decimals=1)\n\n\n\n\n\n\n\n\n\nid\ny\nevent\neventtime\ncenstime\ntrt\nsoc_alt\nhist1\nhist2\nsoc\ntrt_ind\narm\n\n\n\n\n1\n1\n7.7\n0\n12.3\n7.7\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n2\n2\n0.1\n0\n20.4\n0.1\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n3\n3\n4.7\n0\n10.3\n4.7\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n4\n4\n2.7\n0\n19.0\n2.7\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n5\n5\n3.6\n1\n3.6\n7.5\ncontrol\n0\n0\n0\nChemoA\n0\ncontrolChemoA\n\n\n6..199\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n200\n200\n11.8\n0\n17.1\n11.8\nactive\n0\n0\n0\nChemoA\n1\nactiveChemoA\n\n\n\n\n\n\n\ntable(sim$arm)\n\n\n activeChemoA controlChemoA  activeChemoB controlChemoB \n           72            73            28            27 \n\ntable(hdata1$arm)\n\n\n activeChemoA controlChemoA  activeChemoB controlChemoB \n            0           400             0             0",
     "crumbs": [
       "Case studies",
       "<span class='chapter-number'>9</span>  <span class='chapter-title'>Time-to-event data</span>"
@@ -725,7 +725,7 @@
     "href": "src/02i_time_to_event.html#implementationresults",
     "title": "9  Time-to-event data",
     "section": "9.6 Implementation/Results",
-    "text": "9.6 Implementation/Results\n\n9.6.1 Construct customized coded contrasts\nFirst we define customized contrasts which define the parametrization of the model as defined in the Model sepcification section.\nThe mapping of data strata to parameters is:\n\ncc_inv &lt;- matrix(c(1/4,  1/4,  1/4,  1/4,\n                   1/2, -1/2,  1/2, -1/2,\n                   1/2, -1/2, -1/2,  1/2,\n                   0,   -1,    0,    1),\nnrow=4, ncol=4, byrow=TRUE,\ndimnames=list(contrast=c(\"intercept\", \"deltaEffectAvg\",\n                         \"deltaEffect\", \"deltaControl\"),\n              arm=c(\"activeChemoA\", \"controlChemoA\",\n                    \"activeChemoB\", \"controlChemoB\")))\nMASS::fractions(cc_inv)\n\n                arm\ncontrast         activeChemoA controlChemoA activeChemoB controlChemoB\n  intercept       1/4          1/4           1/4          1/4         \n  deltaEffectAvg  1/2         -1/2           1/2         -1/2         \n  deltaEffect     1/2         -1/2          -1/2          1/2         \n  deltaControl      0           -1             0            1         \n\n\nintercept (\\(\\gamma\\)) is the overall mean while deltaControl (\\(\\delta_{\\mbox{control}}\\)) is the difference between the two control arms.\nBy way of defining deltaEffectAvg (\\(\\delta_{\\mbox{EffectAvg}}\\)) we obtain the average treatment effect. Setting deltaEffect (\\(\\delta_{\\mbox{effect}}\\)) to half of the difference between the two treatment effects one arrives at the following symmetric relations to obtain the chemo therapy specific treatment effect:\n\\[\\delta_{\\mbox{EffectAvg}} = \\frac{1}{2} \\, [\n(\\lambda_{\\mbox{activeChemoA}} - \\lambda_{\\mbox{controlChemoA}}) +\n(\\lambda_{\\mbox{activeChemoB}} - \\lambda_{\\mbox{controlChemoB}})]\\]\n\\[\\delta_{\\mbox{effect}} = \\frac{1}{2} \\, [ (\\lambda_{\\mbox{activeChemoA}} - \\lambda_{\\mbox{controlChemoA}}) - (\\lambda_{\\mbox{activeChemoB}} - \\lambda_{\\mbox{controlChemoB}})]\\]\n\\[ \\Leftrightarrow \\delta_{\\mbox{effectChemoA}}  =\n\\delta_{\\mbox{effectAvg}} + \\delta_{\\mbox{effect}}\\]\n\\[ \\Leftrightarrow \\delta_{\\mbox{effectChemoB}}  = \\delta_{\\mbox{effectAvg}} - \\delta_{\\mbox{effect}}\\]\nR requires to define the inverse relationship, i.e. the mapping from parameters to data strata (as needed to define the design matrix). Therefore, the contrast matrix \\(C\\) is:\n\ncc &lt;- solve(cc_inv)\nMASS::fractions(cc)\n\n              intercept deltaEffectAvg deltaEffect deltaControl\nactiveChemoA     1       1/2              1        -1/2        \ncontrolChemoA    1      -1/2              0        -1/2        \nactiveChemoB     1       1/2             -1         1/2        \ncontrolChemoB    1      -1/2              0         1/2        \n\n\nIn order to setup custom contrasts of discrete explanatory variables like treatment arms, R requires to define the explanatory variable as factor, which has been defined accordingly in the simulation routine of this case study:\n\nlevels(sim$arm)\n\n[1] \"activeChemoA\"  \"controlChemoA\" \"activeChemoB\"  \"controlChemoB\"\n\n\nAs R sets up the overall intercept separatley, we need to drop the overall intercept definition from \\(C\\) and define the contrast of the factor variable to be equal to the matrix \\(C\\):\n\ncontrasts(sim$arm) &lt;- cc[,-1]\n\nThen to check if the contrasts work as expected, we can use a dummy fit to see if things get used correctly and the results here are consistent with our definitions:\n\nlfit1 &lt;- lm(y ~ 1 + arm, sim)\nsummary(lfit1)\n\n\nCall:\nlm(formula = y ~ 1 + arm, data = sim)\n\nResiduals:\n   Min     1Q Median     3Q    Max \n-5.506 -2.909 -1.461  1.261 17.064 \n\nCoefficients:\n                  Estimate Std. Error t value Pr(&gt;|t|)    \n(Intercept)         4.5023     0.3349  13.444   &lt;2e-16 ***\narmdeltaEffectAvg   0.9641     0.6698   1.439    0.152    \narmdeltaEffect     -0.4516     0.6698  -0.674    0.501    \narmdeltaControl     0.6574     0.9526   0.690    0.491    \n---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 4.229 on 196 degrees of freedom\nMultiple R-squared:  0.02433,   Adjusted R-squared:  0.009398 \nF-statistic: 1.629 on 3 and 196 DF,  p-value: 0.1839\n\n\nWe can see that above model specification for lm estimates the parameter coefficient for the three contrasts armdeltaEffectAvg (\\(\\delta_{\\mbox{EffectAvg}}\\)), armdeltaEffect (\\(\\delta_{\\mbox{effect}}\\)) and armdeltaControl (\\(\\delta_{\\mbox{control}}\\)).\n\n\n9.6.2 Model fit with brms\nNow we start with analyzing the simulated randomized study and specify the brms model using a Weibull distribution. In a first version, we do not include the literature historical data which requires the definition of a custom design matrix as will be explained below.\nFirst we define with a call to bf the model formula and model family. This defines the model’s likelihood and the parametrization.\n\nmodel_weibull &lt;- bf(y | cens(1-event) ~ 1 + arm, family=weibull())\n\nThe left hand side of the formula, y | cens(1-event) ~ ..., denotes with \\(y\\) the data being modeled - the time variable -, while cens is a response modifier added with a vertical bar| to the left hand side. Here, cens defines which data rows correspond to (right) censored observations. The censoring information is passed into cens by the variable 1-event. The right hand side defines the linear regressor formula used for the mean parameter of the model, which is in this case modeled by default with a \\(\\log\\) link function. The formula 1 + arm defines an overall intercept and the covariate arm (coded as factor in the data) is used with it’s custom constrasts as explanatory variable. Finally, the family argument specifies that the event and censoring times stored in \\(y\\) are distributed according to a weibull probability density distribution.\nWith the model (and data) being defined, we are left to specify the priors for the model. With the help of the get_prior function we can report which default priors and parameters were instantiated by brms and hence need an assignment of a prior.\n\ngt(get_prior(model_weibull, sim))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nHere we define the priors with standard deviations stored in variables, which are passed as stanvars to the brms model. Doing so allows to change the prior quickly, since the Stan model as generated by brms remains unchanged whenever the standard deviation of a prior is changed. A priori it is expected that the outcome for the two control groups will be similar, so the prior on deltaControl (\\(\\delta_{\\mbox{control}}\\)) is centered on \\(0\\). The same is a priori assumed for the deltaEffect (\\(\\delta_{\\mbox{effect}}\\)) parameter.\n\nmodel_weibull_prior &lt;- prior(normal(meanInter, log(4)/1.64),\n                             class=\"Intercept\") +\n  prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +\n  prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +\n  prior(normal(0, sdDeltaControl), coef=armdeltaControl) +\n  prior(gamma(3, 2.7), class=shape)\n\nThe prior for the intercept, \\(\\gamma\\), is set to a mean approximately corresponding to an expected median PFS on anti-PD-1 + SoC, as estimated from CheckMate 649 (here 8 month)\n\nprior_mean &lt;- stanvar(log(8/log(2)), \"meanInter\")\n\nThe standard deviations are setup with reference to a \\(90 \\%\\) credible interval. Hence, \\(log(2)/1.64\\) for the standard deviation of the prior for deltaEffectAvg (\\(\\delta_{\\mbox{effectAvg}}\\)) means that the event process can increase/decrease in progression speed by a factor of \\(2\\) at most.\n\nprior_sd &lt;- stanvar(log(2)/1.64, \"sdEffectAvg\") +\n  stanvar(log(1.25)/1.64, \"sdDeltaEffect\")+\n  stanvar(log(1.25)/1.64, \"sdDeltaControl\")\n\nNow the model and prior specifications are complete and we are ready to run the model in brms. First fit the model using the randomized study data\n\nfit_weibull &lt;- brm(model_weibull, data=sim, prior=model_weibull_prior,\n               stanvars=prior_sd + prior_mean,\n               seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.5 seconds.\nChain 2 finished in 0.6 seconds.\nChain 3 finished in 0.6 seconds.\nChain 4 finished in 0.5 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.5 seconds.\nTotal execution time: 2.7 seconds.\n\n\nWhen looking at the model summary, it is preferable to use robust=TRUE to get the median estimate of the posterior rather than the mean (which would change due to transforms)\n\nsummary(fit_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm \n   Data: sim (Number of observations: 200) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.16      0.12     1.98     2.38 1.00     3998     2848\narmdeltaEffectAvg     0.26      0.19    -0.05     0.58 1.00     3806     2833\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4514     3328\narmdeltaControl       0.06      0.12    -0.13     0.26 1.00     4265     3056\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     0.97      0.08     0.85     1.10 1.00     4014     3098\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSummarize posterior draws based on point estimates, estimation errors and quantiles. And we transform the estimates to the original scale\n\npost_sum &lt;- posterior_summary(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)\nprint(exp(post_sum[,\"Estimate\"]), digits=2)\n\n        b_Intercept b_armdeltaEffectAvg    b_armdeltaEffect   b_armdeltaControl               shape           Intercept \n            8.6e+00             1.3e+00             1.0e+00             1.1e+00             2.6e+00             8.5e+00 \n             lprior                lp__ \n            1.0e+00            5.4e-140 \n\n\nThese estimates are the AFT regression coefficients such that the effect of treatment is to slowdown the time to an event and thus the treatment effect estimate \\(\\delta_{\\mbox{effectAvg}}\\) is positive. To convert to hazard ratios, one needs to multiply with the negative of the estimated shape parameter. Doing so for all 3 regression coefficients by using the rvar facility (provided from the posterior R package) yields:\n\npost_rv &lt;- as_draws_rvars(fit_weibull)\npost_rv_HR &lt;- with(post_rv, exp(-1 * shape * c(b_armdeltaEffectAvg, b_armdeltaEffect, b_armdeltaControl)))\nnames(post_rv_HR) &lt;- c(\"deltaEffectAvg\", \"deltaEffect\", \"deltaControl\")\nprint(summarize_draws(post_rv_HR), digits=2)\n\n# A tibble: 3 × 10\n  variable                    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail\n  &lt;chr&gt;                      &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt;    &lt;dbl&gt;\n1 post_rv_HR[deltaEffectAvg] 0.788  0.777 0.146 0.145 0.568  1.05  1.00    3872.    2876.\n2 post_rv_HR[deltaEffect]    1.00   0.995 0.104 0.103 0.839  1.18  1.00    4515.    3209.\n3 post_rv_HR[deltaControl]   0.950  0.947 0.110 0.107 0.777  1.13  1.00    4282.    2919.\n\n\nObtain posterior mean estimates (mean rate). Median of the posterior estimates indicated a roughly 20% increase in survival in activeChemoA compared to controlChemoA.\n\nnd &lt;- data.frame(y=0, event=0, arm=levels(sim$arm))\npost_mean &lt;- posterior_epred(fit_weibull, newdata=nd)\ncolnames(post_mean) &lt;- nd$arm\nquantile(post_mean[,\"activeChemoA\"] / post_mean[,\"controlChemoA\"],\n         probs=c(0.05, 0.5, 0.95))\n\n       5%       50%       95% \n0.9425351 1.3043160 1.8179356 \n\n\nNow extract posterior coefficients with a 90% credible interval\n\npost_est &lt;- fixef(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)\ngt(as.data.frame(post_est)) %&gt;% fmt_number(everything(), decimals=2)\n\n\n\n\n\n\n\n\nEstimate\nEst.Error\nQ5\nQ95\n\n\n\n\n2.16\n0.12\n1.98\n2.38\n\n\n0.26\n0.19\n−0.05\n0.58\n\n\n0.01\n0.11\n−0.17\n0.18\n\n\n0.06\n0.12\n−0.13\n0.26\n\n\n\n\n\n\n\n\nDefine trial successfull if the deltaEffectAvg exceed 0 (which means the survival time is increased) and the result is not significant here.\n\npost_est[\"armdeltaEffectAvg\",\"Q5\"] &gt; 0\n\n[1] FALSE\n\n\n\n\n9.6.3 Incoporate simulated historical data\nSince we have the simulated historical data hist1 on one control arm, we can incorporate the information by adding the corresponding covariate hist1 (indicator variable being \\(1\\) for the historical study and \\(0\\) otherwise) when specifying the formula for meta-analytic-combined (MAC) analysis:\n\nmodel_mac_weibull &lt;- bf(y | cens(1-event) ~ 1 + arm + hist1, family=weibull())\n\nAnd get_prior function shows that in addition to the priors for the weibull model on randomised data, the co-data model needs one more prior on the historical data.\n\ngt(get_prior(model_mac_weibull, sim))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\n\nb\nhist1\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nWe use a student with 6 degrees of freedom for the regression coefficient for the historical data prior. Using a Student-t distribution robustifies against possible data discordance between the randomised study and the historical data, since very large deviations between the two data sets are admitted by the heavy tails of the prior. The prior is centered at a mean of \\(0\\), corresponding to no expected difference in the outcome for historical and current study.\n\nmodel_weibull_mac_prior &lt;- model_weibull_prior +\n  prior(student_t(6, 0, sdHist), coef=hist1)\n\nA variance allowing substantial differences, a \\(1.8\\) fold of increase/decrease, between the randomised study and historical data is set.\n\nprior_sdHist &lt;- stanvar(log(1.8)/1.64, \"sdHist\")\n\nNow we combine the simulated trial data and historical data to fit the co-data MAC model. Don’t forget to set contrasts of arm in the combined data to the contrast matrix.\n\ncomb_data &lt;- bind_rows(sim, hdata1)\ncontrasts(comb_data$arm) &lt;- cc[,-1]\nfit_mac_weibull &lt;- brm(model_mac_weibull, data=comb_data,\n                       prior=model_weibull_mac_prior,\n                   stanvars=prior_sd + prior_mean + prior_sdHist,\n                   seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.8 seconds.\nChain 2 finished in 1.8 seconds.\nChain 3 finished in 1.7 seconds.\nChain 4 finished in 1.9 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.8 seconds.\nTotal execution time: 7.5 seconds.\n\n\nAnd we can compare the results from two models. Inclusion of historical data on controlChemoA affect the estimates of coefficients for overall mean and EffectAverage.\n\nsummary(fit_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm \n   Data: sim (Number of observations: 200) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.16      0.12     1.98     2.38 1.00     3998     2848\narmdeltaEffectAvg     0.26      0.19    -0.05     0.58 1.00     3806     2833\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4514     3328\narmdeltaControl       0.06      0.12    -0.13     0.26 1.00     4265     3056\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     0.97      0.08     0.85     1.10 1.00     4014     3098\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nsummary(fit_mac_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm + hist1 \n   Data: comb_data (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.30 1.00     4909     3769\narmdeltaEffectAvg     0.29      0.18    -0.01     0.59 1.00     4051     2881\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4748     3261\narmdeltaControl       0.07      0.12    -0.13     0.27 1.00     4343     2769\nhist1                -0.20      0.14    -0.43     0.02 1.00     4039     2944\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.94     1.07 1.00     4712     3074\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nNext, we check the informativeness of the historical data by only including the historical data. Note that we have to request that factor levels for which no realization is present in the data set must be kept, since only the controlChemoA arm is present in the historical data.\n\ncontrasts(hdata1$arm) &lt;- cc[,-1]\nfit_mac_prior_weibull &lt;- brm(model_mac_weibull, data=hdata1,\n                             prior=model_weibull_mac_prior,\n                             drop_unused_levels=FALSE,\n                         stanvars=prior_sd + prior_mean + prior_sdHist,\n                         seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.3 seconds.\nChain 2 finished in 1.1 seconds.\nChain 3 finished in 1.2 seconds.\nChain 4 finished in 1.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.2 seconds.\nTotal execution time: 5.1 seconds.\n\ngt(as.data.frame(fixef(fit_mac_prior_weibull, robust=TRUE))) %&gt;% fmt_number(everything(), decimals=2)\n\n\n\n\n\n\n\n\nEstimate\nEst.Error\nQ2.5\nQ97.5\n\n\n\n\n1.73\n0.47\n0.76\n2.69\n\n\n0.01\n0.42\n−0.80\n0.85\n\n\n0.00\n0.14\n−0.28\n0.27\n\n\n0.00\n0.14\n−0.27\n0.27\n\n\n0.01\n0.38\n−0.86\n0.87\n\n\n\n\n\n\n\n\nWe see that we essentially sample the prior for most regression coefficients, but we do substantially reduce the standard error of the intercept estimate which is almost halved in comparison to the prior standard deviation of 0.85.\n\n\n9.6.4 Including historical data of average SoC\nSo far we have setup a custom parametrization allowing to define priors aligned with prior expectations (small difference of the treatment effect and small difference between SoC). However, the CheckMate 649 historical data reports the effect of an average SoC treatment. Thus, the data of CheckMate 649 is incompatible with the 4 defined categorical factor levels so far as it does not correspond to any of the 4 treatment arms in the current data set. However, as the custom factor contrasts and factor level definitions are translated by R into a design matrix, we can customize the created design matrix directly. This requires to setup the design matrix manually giving us the freedom to define the design matrix appropiatley for the CheckMate 649 data set.\nThere is an alternative way to use the customized contrast and setup analysis by using respective indicator flags, which is created using model.matrix. The model.matrix function is used internally by brms to create the design matrix for the regression problem\n\nXind &lt;- model.matrix(~ 1 + arm, comb_data)\ngt_preview(Xind)\n\n\n\n\n\n\n\n\n\n(Intercept)\narmdeltaEffectAvg\narmdeltaEffect\narmdeltaControl\n\n\n\n\n1\n1\n-0.5\n0\n-0.5\n\n\n2\n1\n0.5\n1\n-0.5\n\n\n3\n1\n-0.5\n0\n-0.5\n\n\n4\n1\n0.5\n1\n-0.5\n\n\n5\n1\n-0.5\n0\n-0.5\n\n\n6..599\n\n\n\n\n\n\n600\n1\n-0.5\n0\n-0.5\n\n\n\n\n\n\n\n\nBy adding the generated design matrix to the original data set, we can then simply delcare the model formula using these created indicator variables:\n\ncomb_data_ind &lt;- cbind(comb_data, Xind[,-1])\n\nThus we use those indicator variables in the fit as if these were continuous covariates which replaces the use of the categorical factor variable arm used in the first approach.\n\nmodel_mac_weibull_ind &lt;- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +\n                            armdeltaEffect + armdeltaControl + hist1,\n                            family=weibull())\n\nAnd query what prior parameters need to be set\n\ngt(get_prior(model_mac_weibull_ind, comb_data_ind))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\n\nb\nhist1\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nFit the model using indicator flags with the same priors as before and compare with the previous one fitted with arm variable, the results are the same.\n\nmodel_weibull_mac_prior_ind &lt;- prior(normal(meanInter, log(4)/1.64),\n                                     class=\"Intercept\") +\n    prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +\n    prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +\n    prior(normal(0, sdDeltaControl), coef=armdeltaControl) +\n    prior(normal(0, sdHist), coef=hist1)\n\nfit_mac_weibull_ind &lt;- brm(model_mac_weibull_ind, data=comb_data_ind,\n                           prior=model_weibull_mac_prior_ind,\n                       stanvars=prior_sd + prior_mean + prior_sdHist,\n                       seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.8 seconds.\nChain 2 finished in 1.7 seconds.\nChain 3 finished in 1.8 seconds.\nChain 4 finished in 1.7 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.8 seconds.\nTotal execution time: 7.3 seconds.\n\nsummary(fit_mac_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm + hist1 \n   Data: comb_data (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.30 1.00     4909     3769\narmdeltaEffectAvg     0.29      0.18    -0.01     0.59 1.00     4051     2881\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4748     3261\narmdeltaControl       0.07      0.12    -0.13     0.27 1.00     4343     2769\nhist1                -0.20      0.14    -0.43     0.02 1.00     4039     2944\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.94     1.07 1.00     4712     3074\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nsummary(fit_mac_weibull_ind, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 \n   Data: comb_data_ind (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.29 1.00     5075     3553\narmdeltaEffectAvg     0.28      0.18     0.00     0.57 1.00     4184     3537\narmdeltaEffect        0.01      0.10    -0.17     0.17 1.00     4688     3177\narmdeltaControl       0.07      0.11    -0.12     0.25 1.00     5187     3251\nhist1                -0.19      0.14    -0.43     0.02 1.00     4334     3137\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.93     1.07 1.00     4680     2946\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nNote that the equality holds exactly for the entire posterior, since the very same problem (using the same random seed) has been fitted by brms:\n\npost_mac &lt;- as_draws_matrix(fit_mac_weibull)\npost_mac_ind &lt;- as_draws_matrix(fit_mac_weibull_ind)\nall(post_mac == post_mac_ind)\n\n[1] FALSE\n\n\nThe benefit of this approach is that we can then include external data corresponding to the controlAverage from the second historical study CheckMate 649, which is not one of the treatment groups but the average effect of control groups, \\(\\frac{1}{2} \\, (\\lambda_{\\mbox{controlChemoA}}+\\lambda_{\\mbox{controlChemoB}})\\).\nTo now find the indicators needed for the effect of controlAverage in the given parametrization so far, we can obtain this via simple algebra as:\n\ncc &lt;- MASS::fractions(solve(cc_inv))\ncc\n\n              intercept deltaEffectAvg deltaEffect deltaControl\nactiveChemoA     1       1/2              1        -1/2        \ncontrolChemoA    1      -1/2              0        -1/2        \nactiveChemoB     1       1/2             -1         1/2        \ncontrolChemoB    1      -1/2              0         1/2        \n\n0.5*(cc[\"controlChemoA\",] + cc[\"controlChemoB\",])\n\n     intercept deltaEffectAvg    deltaEffect   deltaControl \n             1           -1/2              0              0 \n\n\nHence we must set the deltaEffectAvg column to \\(-\\frac{1}{2}\\) and the others to \\(0\\) to get the indicator for controlAverage (the overall intercept is always added to the linear predictor as we define 1 as the first term of the model formula) in order to include the new historical data from CheckMate649.\n\nhdata2_ind &lt;- rename(hdata2, y=time, event=status) %&gt;%\n  mutate(hist1=0, hist2=1, armdeltaEffectAvg=-0.5,\n         armdeltaEffect=0, armdeltaControl=0)\nall_comb_data_ind &lt;- bind_rows(comb_data_ind, hdata2_ind)\n\nTo include the second historical data, we specify the formula and the priors as we did for one historical study.\n\nmodel_mac_weibull_ind_all &lt;- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +\n                              armdeltaEffect + armdeltaControl + hist1+ hist2,\n                              family=weibull())\nmodel_weibull_mac_prior_all &lt;- model_weibull_mac_prior +\n  prior(normal(0, sdHist), coef=hist2)\n\nFit the model and compare:\n\nfit_all_mac_weibull_ind &lt;- brm(model_mac_weibull_ind_all, data=all_comb_data_ind,\n                               prior=model_weibull_mac_prior_all,\n                               stanvars=prior_sd + prior_mean + prior_sdHist,\n                               seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 6.4 seconds.\nChain 2 finished in 6.5 seconds.\nChain 3 finished in 6.4 seconds.\nChain 4 finished in 6.2 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 6.4 seconds.\nTotal execution time: 25.9 seconds.\n\nfit_mac_weibull_ind\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 \n   Data: comb_data_ind (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.93     2.33 1.00     5075     3553\narmdeltaEffectAvg     0.29      0.17    -0.05     0.63 1.00     4184     3537\narmdeltaEffect        0.01      0.10    -0.20     0.21 1.00     4688     3177\narmdeltaControl       0.07      0.11    -0.16     0.29 1.00     5187     3251\nhist1                -0.20      0.14    -0.47     0.06 1.00     4334     3137\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.92     1.09 1.00     4680     2946\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nfit_all_mac_weibull_ind\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 + hist2 \n   Data: all_comb_data_ind (Number of observations: 1389) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.09      0.09     1.92     2.27 1.00     4002     3141\narmdeltaEffectAvg     0.23      0.16    -0.09     0.55 1.00     2743     2920\narmdeltaEffect        0.00      0.10    -0.19     0.20 1.00     3616     2884\narmdeltaControl       0.09      0.12    -0.14     0.32 1.00     3555     3039\nhist1                -0.23      0.13    -0.49     0.00 1.00     2287     2849\nhist2                 0.48      0.12     0.24     0.71 1.00     2108     1868\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.12      0.03     1.07     1.18 1.00     4034     3292\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).",
+    "text": "9.6 Implementation/Results\n\n9.6.1 Construct customized coded contrasts\nFirst we define customized contrasts which define the parametrization of the model as defined in the Model sepcification section.\nThe mapping of data strata to parameters is:\n\ncc_inv &lt;- matrix(c(1/4,  1/4,  1/4,  1/4,\n                   1/2, -1/2,  1/2, -1/2,\n                   1/2, -1/2, -1/2,  1/2,\n                   0,   -1,    0,    1),\nnrow=4, ncol=4, byrow=TRUE,\ndimnames=list(contrast=c(\"intercept\", \"deltaEffectAvg\",\n                         \"deltaEffect\", \"deltaControl\"),\n              arm=c(\"activeChemoA\", \"controlChemoA\",\n                    \"activeChemoB\", \"controlChemoB\")))\nMASS::fractions(cc_inv)\n\n                arm\ncontrast         activeChemoA controlChemoA activeChemoB controlChemoB\n  intercept       1/4          1/4           1/4          1/4         \n  deltaEffectAvg  1/2         -1/2           1/2         -1/2         \n  deltaEffect     1/2         -1/2          -1/2          1/2         \n  deltaControl      0           -1             0            1         \n\n\nintercept (\\(\\gamma\\)) is the overall mean while deltaControl (\\(\\delta_{\\mbox{control}}\\)) is the difference between the two control arms.\nBy way of defining deltaEffectAvg (\\(\\delta_{\\mbox{EffectAvg}}\\)) we obtain the average treatment effect. Setting deltaEffect (\\(\\delta_{\\mbox{effect}}\\)) to half of the difference between the two treatment effects one arrives at the following symmetric relations to obtain the chemo therapy specific treatment effect:\n\\[\\delta_{\\mbox{EffectAvg}} = \\frac{1}{2} \\, [\n(\\lambda_{\\mbox{activeChemoA}} - \\lambda_{\\mbox{controlChemoA}}) +\n(\\lambda_{\\mbox{activeChemoB}} - \\lambda_{\\mbox{controlChemoB}})]\\]\n\\[\\delta_{\\mbox{effect}} = \\frac{1}{2} \\, [ (\\lambda_{\\mbox{activeChemoA}} - \\lambda_{\\mbox{controlChemoA}}) - (\\lambda_{\\mbox{activeChemoB}} - \\lambda_{\\mbox{controlChemoB}})]\\]\n\\[ \\Leftrightarrow \\delta_{\\mbox{effectChemoA}}  =\n\\delta_{\\mbox{effectAvg}} + \\delta_{\\mbox{effect}}\\]\n\\[ \\Leftrightarrow \\delta_{\\mbox{effectChemoB}}  = \\delta_{\\mbox{effectAvg}} - \\delta_{\\mbox{effect}}\\]\nR requires to define the inverse relationship, i.e. the mapping from parameters to data strata (as needed to define the design matrix). Therefore, the contrast matrix \\(C\\) is:\n\ncc &lt;- solve(cc_inv)\nMASS::fractions(cc)\n\n              intercept deltaEffectAvg deltaEffect deltaControl\nactiveChemoA     1       1/2              1        -1/2        \ncontrolChemoA    1      -1/2              0        -1/2        \nactiveChemoB     1       1/2             -1         1/2        \ncontrolChemoB    1      -1/2              0         1/2        \n\n\nIn order to setup custom contrasts of discrete explanatory variables like treatment arms, R requires to define the explanatory variable as factor, which has been defined accordingly in the simulation routine of this case study:\n\nlevels(sim$arm)\n\n[1] \"activeChemoA\"  \"controlChemoA\" \"activeChemoB\"  \"controlChemoB\"\n\n\nAs R sets up the overall intercept separatley, we need to drop the overall intercept definition from \\(C\\) and define the contrast of the factor variable to be equal to the matrix \\(C\\):\n\ncontrasts(sim$arm) &lt;- cc[,-1]\n\nThen to check if the contrasts work as expected, we can use a dummy fit to see if things get used correctly and the results here are consistent with our definitions:\n\nlfit1 &lt;- lm(y ~ 1 + arm, sim)\nsummary(lfit1)\n\n\nCall:\nlm(formula = y ~ 1 + arm, data = sim)\n\nResiduals:\n   Min     1Q Median     3Q    Max \n-5.506 -2.909 -1.461  1.261 17.064 \n\nCoefficients:\n                  Estimate Std. Error t value Pr(&gt;|t|)    \n(Intercept)         4.5023     0.3349  13.444   &lt;2e-16 ***\narmdeltaEffectAvg   0.9641     0.6698   1.439    0.152    \narmdeltaEffect     -0.4516     0.6698  -0.674    0.501    \narmdeltaControl     0.6574     0.9526   0.690    0.491    \n---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 4.229 on 196 degrees of freedom\nMultiple R-squared:  0.02433,   Adjusted R-squared:  0.009398 \nF-statistic: 1.629 on 3 and 196 DF,  p-value: 0.1839\n\n\nWe can see that above model specification for lm estimates the parameter coefficient for the three contrasts armdeltaEffectAvg (\\(\\delta_{\\mbox{EffectAvg}}\\)), armdeltaEffect (\\(\\delta_{\\mbox{effect}}\\)) and armdeltaControl (\\(\\delta_{\\mbox{control}}\\)).\n\n\n9.6.2 Model fit with brms\nNow we start with analyzing the simulated randomized study and specify the brms model using a Weibull distribution. In a first version, we do not include the literature historical data which requires the definition of a custom design matrix as will be explained below.\nFirst we define with a call to bf the model formula and model family. This defines the model’s likelihood and the parametrization.\n\nmodel_weibull &lt;- bf(y | cens(1-event) ~ 1 + arm, family=weibull())\n\nThe left hand side of the formula, y | cens(1-event) ~ ..., denotes with \\(y\\) the data being modeled - the time variable -, while cens is a response modifier added with a vertical bar| to the left hand side. Here, cens defines which data rows correspond to (right) censored observations. The censoring information is passed into cens by the variable 1-event. The right hand side defines the linear regressor formula used for the mean parameter of the model, which is in this case modeled by default with a \\(\\log\\) link function. The formula 1 + arm defines an overall intercept and the covariate arm (coded as factor in the data) is used with it’s custom constrasts as explanatory variable. Finally, the family argument specifies that the event and censoring times stored in \\(y\\) are distributed according to a weibull probability density distribution.\nWith the model (and data) being defined, we are left to specify the priors for the model. With the help of the get_prior function we can report which default priors and parameters were instantiated by brms and hence need an assignment of a prior.\n\ngt(get_prior(model_weibull, sim))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nHere we define the priors with standard deviations stored in variables, which are passed as stanvars to the brms model. Doing so allows to change the prior quickly, since the Stan model as generated by brms remains unchanged whenever the standard deviation of a prior is changed. A priori it is expected that the outcome for the two control groups will be similar, so the prior on deltaControl (\\(\\delta_{\\mbox{control}}\\)) is centered on \\(0\\). The same is a priori assumed for the deltaEffect (\\(\\delta_{\\mbox{effect}}\\)) parameter.\n\nmodel_weibull_prior &lt;- prior(normal(meanInter, log(4)/1.64),\n                             class=\"Intercept\") +\n  prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +\n  prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +\n  prior(normal(0, sdDeltaControl), coef=armdeltaControl) +\n  prior(gamma(3, 2.7), class=shape)\n\nThe prior for the intercept, \\(\\gamma\\), is set to a mean approximately corresponding to an expected median PFS on anti-PD-1 + SoC, as estimated from CheckMate 649 (here 8 month)\n\nprior_mean &lt;- stanvar(log(8/log(2)), \"meanInter\")\n\nThe standard deviations are setup with reference to a \\(90 \\%\\) credible interval. Hence, \\(log(2)/1.64\\) for the standard deviation of the prior for deltaEffectAvg (\\(\\delta_{\\mbox{effectAvg}}\\)) means that the event process can increase/decrease in progression speed by a factor of \\(2\\) at most.\n\nprior_sd &lt;- stanvar(log(2)/1.64, \"sdEffectAvg\") +\n  stanvar(log(1.25)/1.64, \"sdDeltaEffect\")+\n  stanvar(log(1.25)/1.64, \"sdDeltaControl\")\n\nNow the model and prior specifications are complete and we are ready to run the model in brms. First fit the model using the randomized study data\n\nfit_weibull &lt;- brm(model_weibull, data=sim, prior=model_weibull_prior,\n               stanvars=prior_sd + prior_mean,\n               seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.5 seconds.\nChain 2 finished in 0.6 seconds.\nChain 3 finished in 0.5 seconds.\nChain 4 finished in 0.5 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.5 seconds.\nTotal execution time: 2.6 seconds.\n\n\nWhen looking at the model summary, it is preferable to use robust=TRUE to get the median estimate of the posterior rather than the mean (which would change due to transforms)\n\nsummary(fit_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm \n   Data: sim (Number of observations: 200) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.16      0.12     1.98     2.38 1.00     3998     2848\narmdeltaEffectAvg     0.26      0.19    -0.05     0.58 1.00     3806     2833\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4514     3328\narmdeltaControl       0.06      0.12    -0.13     0.26 1.00     4265     3056\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     0.97      0.08     0.85     1.10 1.00     4014     3098\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nSummarize posterior draws based on point estimates, estimation errors and quantiles. And we transform the estimates to the original scale\n\npost_sum &lt;- posterior_summary(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)\nprint(exp(post_sum[,\"Estimate\"]), digits=2)\n\n        b_Intercept b_armdeltaEffectAvg    b_armdeltaEffect   b_armdeltaControl               shape           Intercept \n            8.6e+00             1.3e+00             1.0e+00             1.1e+00             2.6e+00             8.5e+00 \n             lprior                lp__ \n            1.0e+00            5.4e-140 \n\n\nThese estimates are the AFT regression coefficients such that the effect of treatment is to slowdown the time to an event and thus the treatment effect estimate \\(\\delta_{\\mbox{effectAvg}}\\) is positive. To convert to hazard ratios, one needs to multiply with the negative of the estimated shape parameter. Doing so for all 3 regression coefficients by using the rvar facility (provided from the posterior R package) yields:\n\npost_rv &lt;- as_draws_rvars(fit_weibull)\npost_rv_HR &lt;- with(post_rv, exp(-1 * shape * c(b_armdeltaEffectAvg, b_armdeltaEffect, b_armdeltaControl)))\nnames(post_rv_HR) &lt;- c(\"deltaEffectAvg\", \"deltaEffect\", \"deltaControl\")\nprint(summarize_draws(post_rv_HR), digits=2)\n\n# A tibble: 3 × 10\n  variable                    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail\n  &lt;chr&gt;                      &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt;    &lt;dbl&gt;\n1 post_rv_HR[deltaEffectAvg] 0.788  0.777 0.146 0.145 0.568  1.05  1.00    3872.    2876.\n2 post_rv_HR[deltaEffect]    1.00   0.995 0.104 0.103 0.839  1.18  1.00    4515.    3209.\n3 post_rv_HR[deltaControl]   0.950  0.947 0.110 0.107 0.777  1.13  1.00    4282.    2919.\n\n\nObtain posterior mean estimates (mean rate). Median of the posterior estimates indicated a roughly 20% increase in survival in activeChemoA compared to controlChemoA.\n\nnd &lt;- data.frame(y=0, event=0, arm=levels(sim$arm))\npost_mean &lt;- posterior_epred(fit_weibull, newdata=nd)\ncolnames(post_mean) &lt;- nd$arm\nquantile(post_mean[,\"activeChemoA\"] / post_mean[,\"controlChemoA\"],\n         probs=c(0.05, 0.5, 0.95))\n\n       5%       50%       95% \n0.9425351 1.3043160 1.8179356 \n\n\nNow extract posterior coefficients with a 90% credible interval\n\npost_est &lt;- fixef(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)\ngt(as.data.frame(post_est)) %&gt;% fmt_number(everything(), decimals=2)\n\n\n\n\n\n\n\n\nEstimate\nEst.Error\nQ5\nQ95\n\n\n\n\n2.16\n0.12\n1.98\n2.38\n\n\n0.26\n0.19\n−0.05\n0.58\n\n\n0.01\n0.11\n−0.17\n0.18\n\n\n0.06\n0.12\n−0.13\n0.26\n\n\n\n\n\n\n\n\nDefine trial successfull if the deltaEffectAvg exceed 0 (which means the survival time is increased) and the result is not significant here.\n\npost_est[\"armdeltaEffectAvg\",\"Q5\"] &gt; 0\n\n[1] FALSE\n\n\n\n\n9.6.3 Incoporate simulated historical data\nSince we have the simulated historical data hist1 on one control arm, we can incorporate the information by adding the corresponding covariate hist1 (indicator variable being \\(1\\) for the historical study and \\(0\\) otherwise) when specifying the formula for meta-analytic-combined (MAC) analysis:\n\nmodel_mac_weibull &lt;- bf(y | cens(1-event) ~ 1 + arm + hist1, family=weibull())\n\nAnd get_prior function shows that in addition to the priors for the weibull model on randomised data, the co-data model needs one more prior on the historical data.\n\ngt(get_prior(model_mac_weibull, sim))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\n\nb\nhist1\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nWe use a student with 6 degrees of freedom for the regression coefficient for the historical data prior. Using a Student-t distribution robustifies against possible data discordance between the randomised study and the historical data, since very large deviations between the two data sets are admitted by the heavy tails of the prior. The prior is centered at a mean of \\(0\\), corresponding to no expected difference in the outcome for historical and current study.\n\nmodel_weibull_mac_prior &lt;- model_weibull_prior +\n  prior(student_t(6, 0, sdHist), coef=hist1)\n\nA variance allowing substantial differences, a \\(1.8\\) fold of increase/decrease, between the randomised study and historical data is set.\n\nprior_sdHist &lt;- stanvar(log(1.8)/1.64, \"sdHist\")\n\nNow we combine the simulated trial data and historical data to fit the co-data MAC model. Don’t forget to set contrasts of arm in the combined data to the contrast matrix.\n\ncomb_data &lt;- bind_rows(sim, hdata1)\ncontrasts(comb_data$arm) &lt;- cc[,-1]\nfit_mac_weibull &lt;- brm(model_mac_weibull, data=comb_data,\n                       prior=model_weibull_mac_prior,\n                   stanvars=prior_sd + prior_mean + prior_sdHist,\n                   seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.7 seconds.\nChain 2 finished in 1.8 seconds.\nChain 3 finished in 1.7 seconds.\nChain 4 finished in 1.8 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.8 seconds.\nTotal execution time: 7.4 seconds.\n\n\nAnd we can compare the results from two models. Inclusion of historical data on controlChemoA affect the estimates of coefficients for overall mean and EffectAverage.\n\nsummary(fit_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm \n   Data: sim (Number of observations: 200) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.16      0.12     1.98     2.38 1.00     3998     2848\narmdeltaEffectAvg     0.26      0.19    -0.05     0.58 1.00     3806     2833\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4514     3328\narmdeltaControl       0.06      0.12    -0.13     0.26 1.00     4265     3056\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     0.97      0.08     0.85     1.10 1.00     4014     3098\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nsummary(fit_mac_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm + hist1 \n   Data: comb_data (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.30 1.00     4909     3769\narmdeltaEffectAvg     0.29      0.18    -0.01     0.59 1.00     4051     2881\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4748     3261\narmdeltaControl       0.07      0.12    -0.13     0.27 1.00     4343     2769\nhist1                -0.20      0.14    -0.43     0.02 1.00     4039     2944\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.94     1.07 1.00     4712     3074\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nNext, we check the informativeness of the historical data by only including the historical data. Note that we have to request that factor levels for which no realization is present in the data set must be kept, since only the controlChemoA arm is present in the historical data.\n\ncontrasts(hdata1$arm) &lt;- cc[,-1]\nfit_mac_prior_weibull &lt;- brm(model_mac_weibull, data=hdata1,\n                             prior=model_weibull_mac_prior,\n                             drop_unused_levels=FALSE,\n                         stanvars=prior_sd + prior_mean + prior_sdHist,\n                         seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.2 seconds.\nChain 2 finished in 1.1 seconds.\nChain 3 finished in 1.2 seconds.\nChain 4 finished in 1.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.1 seconds.\nTotal execution time: 4.9 seconds.\n\ngt(as.data.frame(fixef(fit_mac_prior_weibull, robust=TRUE))) %&gt;% fmt_number(everything(), decimals=2)\n\n\n\n\n\n\n\n\nEstimate\nEst.Error\nQ2.5\nQ97.5\n\n\n\n\n1.73\n0.47\n0.76\n2.69\n\n\n0.01\n0.42\n−0.80\n0.85\n\n\n0.00\n0.14\n−0.28\n0.27\n\n\n0.00\n0.14\n−0.27\n0.27\n\n\n0.01\n0.38\n−0.86\n0.87\n\n\n\n\n\n\n\n\nWe see that we essentially sample the prior for most regression coefficients, but we do substantially reduce the standard error of the intercept estimate which is almost halved in comparison to the prior standard deviation of 0.85.\n\n\n9.6.4 Including historical data of average SoC\nSo far we have setup a custom parametrization allowing to define priors aligned with prior expectations (small difference of the treatment effect and small difference between SoC). However, the CheckMate 649 historical data reports the effect of an average SoC treatment. Thus, the data of CheckMate 649 is incompatible with the 4 defined categorical factor levels so far as it does not correspond to any of the 4 treatment arms in the current data set. However, as the custom factor contrasts and factor level definitions are translated by R into a design matrix, we can customize the created design matrix directly. This requires to setup the design matrix manually giving us the freedom to define the design matrix appropiatley for the CheckMate 649 data set.\nThere is an alternative way to use the customized contrast and setup analysis by using respective indicator flags, which is created using model.matrix. The model.matrix function is used internally by brms to create the design matrix for the regression problem\n\nXind &lt;- model.matrix(~ 1 + arm, comb_data)\ngt_preview(Xind)\n\n\n\n\n\n\n\n\n\n(Intercept)\narmdeltaEffectAvg\narmdeltaEffect\narmdeltaControl\n\n\n\n\n1\n1\n-0.5\n0\n-0.5\n\n\n2\n1\n0.5\n1\n-0.5\n\n\n3\n1\n-0.5\n0\n-0.5\n\n\n4\n1\n0.5\n1\n-0.5\n\n\n5\n1\n-0.5\n0\n-0.5\n\n\n6..599\n\n\n\n\n\n\n600\n1\n-0.5\n0\n-0.5\n\n\n\n\n\n\n\n\nBy adding the generated design matrix to the original data set, we can then simply delcare the model formula using these created indicator variables:\n\ncomb_data_ind &lt;- cbind(comb_data, Xind[,-1])\n\nThus we use those indicator variables in the fit as if these were continuous covariates which replaces the use of the categorical factor variable arm used in the first approach.\n\nmodel_mac_weibull_ind &lt;- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +\n                            armdeltaEffect + armdeltaControl + hist1,\n                            family=weibull())\n\nAnd query what prior parameters need to be set\n\ngt(get_prior(model_mac_weibull_ind, comb_data_ind))\n\n\n\n\n\n\n\n\nprior\nclass\ncoef\ngroup\nresp\ndpar\nnlpar\nlb\nub\nsource\n\n\n\n\n\nb\n\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaControl\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffect\n\n\n\n\n\n\ndefault\n\n\n\nb\narmdeltaEffectAvg\n\n\n\n\n\n\ndefault\n\n\n\nb\nhist1\n\n\n\n\n\n\ndefault\n\n\nstudent_t(3, 1, 2.5)\nIntercept\n\n\n\n\n\n\n\ndefault\n\n\ngamma(0.01, 0.01)\nshape\n\n\n\n\n\n0\n\ndefault\n\n\n\n\n\n\n\n\nFit the model using indicator flags with the same priors as before and compare with the previous one fitted with arm variable, the results are the same.\n\nmodel_weibull_mac_prior_ind &lt;- prior(normal(meanInter, log(4)/1.64),\n                                     class=\"Intercept\") +\n    prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +\n    prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +\n    prior(normal(0, sdDeltaControl), coef=armdeltaControl) +\n    prior(normal(0, sdHist), coef=hist1)\n\nfit_mac_weibull_ind &lt;- brm(model_mac_weibull_ind, data=comb_data_ind,\n                           prior=model_weibull_mac_prior_ind,\n                       stanvars=prior_sd + prior_mean + prior_sdHist,\n                       seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 1.8 seconds.\nChain 2 finished in 1.7 seconds.\nChain 3 finished in 1.8 seconds.\nChain 4 finished in 1.7 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.7 seconds.\nTotal execution time: 7.2 seconds.\n\nsummary(fit_mac_weibull, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + arm + hist1 \n   Data: comb_data (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.30 1.00     4909     3769\narmdeltaEffectAvg     0.29      0.18    -0.01     0.59 1.00     4051     2881\narmdeltaEffect        0.01      0.11    -0.17     0.18 1.00     4748     3261\narmdeltaControl       0.07      0.12    -0.13     0.27 1.00     4343     2769\nhist1                -0.20      0.14    -0.43     0.02 1.00     4039     2944\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.94     1.07 1.00     4712     3074\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nsummary(fit_mac_weibull_ind, prob=0.9, robust=TRUE)\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 \n   Data: comb_data_ind (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.96     2.29 1.00     5075     3553\narmdeltaEffectAvg     0.28      0.18     0.00     0.57 1.00     4184     3537\narmdeltaEffect        0.01      0.10    -0.17     0.17 1.00     4688     3177\narmdeltaControl       0.07      0.11    -0.12     0.25 1.00     5187     3251\nhist1                -0.19      0.14    -0.43     0.02 1.00     4334     3137\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.93     1.07 1.00     4680     2946\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\n\nNote that the equality holds exactly for the entire posterior, since the very same problem (using the same random seed) has been fitted by brms:\n\npost_mac &lt;- as_draws_matrix(fit_mac_weibull)\npost_mac_ind &lt;- as_draws_matrix(fit_mac_weibull_ind)\nall(post_mac == post_mac_ind)\n\n[1] FALSE\n\n\nThe benefit of this approach is that we can then include external data corresponding to the controlAverage from the second historical study CheckMate 649, which is not one of the treatment groups but the average effect of control groups, \\(\\frac{1}{2} \\, (\\lambda_{\\mbox{controlChemoA}}+\\lambda_{\\mbox{controlChemoB}})\\).\nTo now find the indicators needed for the effect of controlAverage in the given parametrization so far, we can obtain this via simple algebra as:\n\ncc &lt;- MASS::fractions(solve(cc_inv))\ncc\n\n              intercept deltaEffectAvg deltaEffect deltaControl\nactiveChemoA     1       1/2              1        -1/2        \ncontrolChemoA    1      -1/2              0        -1/2        \nactiveChemoB     1       1/2             -1         1/2        \ncontrolChemoB    1      -1/2              0         1/2        \n\n0.5*(cc[\"controlChemoA\",] + cc[\"controlChemoB\",])\n\n     intercept deltaEffectAvg    deltaEffect   deltaControl \n             1           -1/2              0              0 \n\n\nHence we must set the deltaEffectAvg column to \\(-\\frac{1}{2}\\) and the others to \\(0\\) to get the indicator for controlAverage (the overall intercept is always added to the linear predictor as we define 1 as the first term of the model formula) in order to include the new historical data from CheckMate649.\n\nhdata2_ind &lt;- rename(hdata2, y=time, event=status) %&gt;%\n  mutate(hist1=0, hist2=1, armdeltaEffectAvg=-0.5,\n         armdeltaEffect=0, armdeltaControl=0)\nall_comb_data_ind &lt;- bind_rows(comb_data_ind, hdata2_ind)\n\nTo include the second historical data, we specify the formula and the priors as we did for one historical study.\n\nmodel_mac_weibull_ind_all &lt;- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +\n                              armdeltaEffect + armdeltaControl + hist1+ hist2,\n                              family=weibull())\nmodel_weibull_mac_prior_all &lt;- model_weibull_mac_prior +\n  prior(normal(0, sdHist), coef=hist2)\n\nFit the model and compare:\n\nfit_all_mac_weibull_ind &lt;- brm(model_mac_weibull_ind_all, data=all_comb_data_ind,\n                               prior=model_weibull_mac_prior_all,\n                               stanvars=prior_sd + prior_mean + prior_sdHist,\n                               seed=43534, refresh=0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 6.3 seconds.\nChain 2 finished in 6.4 seconds.\nChain 3 finished in 6.3 seconds.\nChain 4 finished in 6.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 6.3 seconds.\nTotal execution time: 25.4 seconds.\n\nfit_mac_weibull_ind\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 \n   Data: comb_data_ind (Number of observations: 600) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.12      0.10     1.93     2.33 1.00     5075     3553\narmdeltaEffectAvg     0.29      0.17    -0.05     0.63 1.00     4184     3537\narmdeltaEffect        0.01      0.10    -0.20     0.21 1.00     4688     3177\narmdeltaControl       0.07      0.11    -0.16     0.29 1.00     5187     3251\nhist1                -0.20      0.14    -0.47     0.06 1.00     4334     3137\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.00      0.04     0.92     1.09 1.00     4680     2946\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).\n\nfit_all_mac_weibull_ind\n\n Family: weibull \n  Links: mu = log; shape = identity \nFormula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 + hist2 \n   Data: all_comb_data_ind (Number of observations: 1389) \n  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;\n         total post-warmup draws = 4000\n\nRegression Coefficients:\n                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nIntercept             2.09      0.09     1.92     2.27 1.00     4002     3141\narmdeltaEffectAvg     0.23      0.16    -0.09     0.55 1.00     2743     2920\narmdeltaEffect        0.00      0.10    -0.19     0.20 1.00     3616     2884\narmdeltaControl       0.09      0.12    -0.14     0.32 1.00     3555     3039\nhist1                -0.23      0.13    -0.49     0.00 1.00     2287     2849\nhist2                 0.48      0.12     0.24     0.71 1.00     2108     1868\n\nFurther Distributional Parameters:\n      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS\nshape     1.12      0.03     1.07     1.18 1.00     4034     3292\n\nDraws were sampled using sample(hmc). For each parameter, Bulk_ESS\nand Tail_ESS are effective sample size measures, and Rhat is the potential\nscale reduction factor on split chains (at convergence, Rhat = 1).",
     "crumbs": [
       "Case studies",
       "<span class='chapter-number'>9</span>  <span class='chapter-title'>Time-to-event data</span>"
diff --git a/sitemap.xml b/sitemap.xml
index ed41a6e..df1e762 100644
--- a/sitemap.xml
+++ b/sitemap.xml
@@ -2,58 +2,58 @@
 <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9">
   <url>
     <loc>https://opensource.nibr.com/RBesT/index.html</loc>
-    <lastmod>2024-04-18T18:57:00.885Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.399Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/01a_introduction.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/01b_basic_workflow.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02_case_studies.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02a_meta_analysis.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02ab_meta_analysis_trtdiff.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02b_dose_finding.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02c_dose_escalation.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02e_multiple_imputation.html</loc>
-    <lastmod>2024-04-18T18:57:00.893Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.403Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02h_mmrm.html</loc>
-    <lastmod>2024-04-18T18:57:00.897Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.407Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02i_time_to_event.html</loc>
-    <lastmod>2024-04-18T18:57:00.897Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.407Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/02j_network_meta_analysis.html</loc>
-    <lastmod>2024-04-18T18:57:00.897Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.407Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/03a_stan_code.html</loc>
-    <lastmod>2024-04-18T18:57:00.897Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.407Z</lastmod>
   </url>
   <url>
     <loc>https://opensource.nibr.com/RBesT/src/03b_parallel.html</loc>
-    <lastmod>2024-04-18T18:57:00.897Z</lastmod>
+    <lastmod>2024-04-18T22:31:52.407Z</lastmod>
   </url>
 </urlset>
diff --git a/src/01b_basic_workflow.html b/src/01b_basic_workflow.html
index 54da709..4eddaf2 100644
--- a/src/01b_basic_workflow.html
+++ b/src/01b_basic_workflow.html
@@ -680,7 +680,7 @@ <h2 data-number="2.5" class="anchored" data-anchor-id="fitting-the-model"><span
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 0.2 seconds.
+Chain 1 finished in 0.1 seconds.
 Chain 2 finished in 0.1 seconds.
 Chain 3 finished in 0.1 seconds.
 Chain 4 finished in 0.1 seconds.
@@ -821,9 +821,9 @@ <h2 data-number="2.6" class="anchored" data-anchor-id="model-summary"><span clas
 <pre><code>Running MCMC with 4 sequential chains...
 
 Chain 1 finished in 0.0 seconds.
-Chain 2 finished in 0.1 seconds.
+Chain 2 finished in 0.0 seconds.
 Chain 3 finished in 0.0 seconds.
-Chain 4 finished in 0.1 seconds.
+Chain 4 finished in 0.0 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 0.0 seconds.
diff --git a/src/02b_dose_finding.html b/src/02b_dose_finding.html
index 048d780..a0542f9 100644
--- a/src/02b_dose_finding.html
+++ b/src/02b_dose_finding.html
@@ -852,35 +852,35 @@ <h3 data-number="5.6.3" class="anchored" data-anchor-id="perform-bayesian-model-
 <pre><code>Running MCMC with 4 sequential chains...
 
 Chain 1 finished in 0.6 seconds.
-Chain 2 finished in 0.8 seconds.
+Chain 2 finished in 0.7 seconds.
 Chain 3 finished in 0.8 seconds.
 Chain 4 finished in 0.8 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 0.7 seconds.
-Total execution time: 3.3 seconds.
+Total execution time: 3.2 seconds.
 
 Running MCMC with 4 sequential chains...
 
 Chain 1 finished in 0.3 seconds.
 Chain 2 finished in 0.2 seconds.
 Chain 3 finished in 0.3 seconds.
-Chain 4 finished in 0.5 seconds.
+Chain 4 finished in 0.4 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 0.3 seconds.
-Total execution time: 1.8 seconds.
+Total execution time: 1.5 seconds.
 
 Running MCMC with 4 sequential chains...
 
 Chain 1 finished in 0.2 seconds.
 Chain 2 finished in 0.4 seconds.
-Chain 3 finished in 0.4 seconds.
+Chain 3 finished in 0.3 seconds.
 Chain 4 finished in 0.3 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 0.3 seconds.
-Total execution time: 1.8 seconds.
+Total execution time: 1.6 seconds.
 
 Running MCMC with 4 sequential chains...
 
@@ -917,7 +917,7 @@ <h3 data-number="5.6.3" class="anchored" data-anchor-id="perform-bayesian-model-
 
 Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 0.3 seconds.
+Chain 1 finished in 0.2 seconds.
 Chain 2 finished in 0.3 seconds.
 Chain 3 finished in 0.3 seconds.
 Chain 4 finished in 0.2 seconds.
diff --git a/src/02e_multiple_imputation.html b/src/02e_multiple_imputation.html
index 0f5bab4..ec42b74 100644
--- a/src/02e_multiple_imputation.html
+++ b/src/02e_multiple_imputation.html
@@ -616,14 +616,14 @@ <h3 data-number="7.5.1" class="anchored" data-anchor-id="imputation-from-observe
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 11.3 seconds.
+Chain 1 finished in 11.2 seconds.
 Chain 2 finished in 10.6 seconds.
-Chain 3 finished in 13.7 seconds.
-Chain 4 finished in 13.2 seconds.
+Chain 3 finished in 13.5 seconds.
+Chain 4 finished in 13.1 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 12.2 seconds.
-Total execution time: 49.2 seconds.</code></pre>
+Mean chain execution time: 12.1 seconds.
+Total execution time: 48.7 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb6"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a><span class="co"># Predictions for missing (lost-to-follow-up) post-treatment values</span></span>
 <span id="cb6-2"><a href="#cb6-2" aria-hidden="true" tabindex="-1"></a><span class="co"># Output format is a matrix with each column being a patient record and each row being one imputation</span></span>
@@ -667,14 +667,14 @@ <h3 data-number="7.5.2" class="anchored" data-anchor-id="imputation-based-on-the
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 10.3 seconds.
-Chain 2 finished in 10.2 seconds.
-Chain 3 finished in 12.9 seconds.
-Chain 4 finished in 12.3 seconds.
+Chain 1 finished in 10.1 seconds.
+Chain 2 finished in 10.1 seconds.
+Chain 3 finished in 12.8 seconds.
+Chain 4 finished in 12.1 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 11.4 seconds.
-Total execution time: 46.3 seconds.</code></pre>
+Mean chain execution time: 11.3 seconds.
+Total execution time: 45.5 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb10"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a>imputations2 <span class="ot">&lt;-</span> <span class="fu">predict</span>(brmfit2, </span>
 <span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a>                        <span class="at">newdata=</span>count_data <span class="sc">%&gt;%</span> </span>
@@ -712,14 +712,14 @@ <h3 data-number="7.5.3" class="anchored" data-anchor-id="imputation-under-a-hypo
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 12.6 seconds.
-Chain 2 finished in 9.7 seconds.
-Chain 3 finished in 10.4 seconds.
-Chain 4 finished in 13.0 seconds.
+Chain 1 finished in 12.4 seconds.
+Chain 2 finished in 9.6 seconds.
+Chain 3 finished in 10.2 seconds.
+Chain 4 finished in 12.8 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 11.4 seconds.
-Total execution time: 46.0 seconds.</code></pre>
+Mean chain execution time: 11.3 seconds.
+Total execution time: 45.4 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb14"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb14-1"><a href="#cb14-1" aria-hidden="true" tabindex="-1"></a>imputations3 <span class="ot">&lt;-</span> <span class="fu">predict</span>(brmfit3, </span>
 <span id="cb14-2"><a href="#cb14-2" aria-hidden="true" tabindex="-1"></a>                        <span class="at">newdata=</span>count_data <span class="sc">%&gt;%</span> </span>
diff --git a/src/02i_time_to_event.html b/src/02i_time_to_event.html
index 113efe8..45011b9 100644
--- a/src/02i_time_to_event.html
+++ b/src/02i_time_to_event.html
@@ -598,14 +598,14 @@ <h3 data-number="9.4.1" class="anchored" data-anchor-id="historical-data"><span
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 1.5 seconds.
+Chain 1 finished in 1.3 seconds.
 Chain 2 finished in 1.4 seconds.
 Chain 3 finished in 1.4 seconds.
 Chain 4 finished in 1.3 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 1.4 seconds.
-Total execution time: 5.9 seconds.</code></pre>
+Mean chain execution time: 1.3 seconds.
+Total execution time: 5.7 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>fit_hist_checkmate</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
@@ -3512,12 +3512,12 @@ <h3 data-number="9.6.2" class="anchored" data-anchor-id="model-fit-with-brms"><s
 
 Chain 1 finished in 0.5 seconds.
 Chain 2 finished in 0.6 seconds.
-Chain 3 finished in 0.6 seconds.
+Chain 3 finished in 0.5 seconds.
 Chain 4 finished in 0.5 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 0.5 seconds.
-Total execution time: 2.7 seconds.</code></pre>
+Total execution time: 2.6 seconds.</code></pre>
 </div>
 </div>
 <p>When looking at the model summary, it is preferable to use <code>robust=TRUE</code> to get the median estimate of the posterior rather than the mean (which would change due to transforms)</p>
@@ -4660,14 +4660,14 @@ <h3 data-number="9.6.3" class="anchored" data-anchor-id="incoporate-simulated-hi
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 1.8 seconds.
+Chain 1 finished in 1.7 seconds.
 Chain 2 finished in 1.8 seconds.
 Chain 3 finished in 1.7 seconds.
-Chain 4 finished in 1.9 seconds.
+Chain 4 finished in 1.8 seconds.
 
 All 4 chains finished successfully.
 Mean chain execution time: 1.8 seconds.
-Total execution time: 7.5 seconds.</code></pre>
+Total execution time: 7.4 seconds.</code></pre>
 </div>
 </div>
 <p>And we can compare the results from two models. Inclusion of historical data on <code>controlChemoA</code> affect the estimates of coefficients for overall mean and <code>EffectAverage</code>.</p>
@@ -4733,14 +4733,14 @@ <h3 data-number="9.6.3" class="anchored" data-anchor-id="incoporate-simulated-hi
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 1.3 seconds.
+Chain 1 finished in 1.2 seconds.
 Chain 2 finished in 1.1 seconds.
 Chain 3 finished in 1.2 seconds.
 Chain 4 finished in 1.1 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 1.2 seconds.
-Total execution time: 5.1 seconds.</code></pre>
+Mean chain execution time: 1.1 seconds.
+Total execution time: 4.9 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb53"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb53-1"><a href="#cb53-1" aria-hidden="true" tabindex="-1"></a><span class="fu">gt</span>(<span class="fu">as.data.frame</span>(<span class="fu">fixef</span>(fit_mac_prior_weibull, <span class="at">robust=</span><span class="cn">TRUE</span>))) <span class="sc">%&gt;%</span> <span class="fu">fmt_number</span>(<span class="fu">everything</span>(), <span class="at">decimals=</span><span class="dv">2</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output-display">
@@ -6337,8 +6337,8 @@ <h3 data-number="9.6.4" class="anchored" data-anchor-id="including-historical-da
 Chain 4 finished in 1.7 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 1.8 seconds.
-Total execution time: 7.3 seconds.</code></pre>
+Mean chain execution time: 1.7 seconds.
+Total execution time: 7.2 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb60"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb60-1"><a href="#cb60-1" aria-hidden="true" tabindex="-1"></a><span class="fu">summary</span>(fit_mac_weibull, <span class="at">prob=</span><span class="fl">0.9</span>, <span class="at">robust=</span><span class="cn">TRUE</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">
@@ -6442,14 +6442,14 @@ <h3 data-number="9.6.4" class="anchored" data-anchor-id="including-historical-da
 <div class="cell-output cell-output-stdout">
 <pre><code>Running MCMC with 4 sequential chains...
 
-Chain 1 finished in 6.4 seconds.
-Chain 2 finished in 6.5 seconds.
-Chain 3 finished in 6.4 seconds.
-Chain 4 finished in 6.2 seconds.
+Chain 1 finished in 6.3 seconds.
+Chain 2 finished in 6.4 seconds.
+Chain 3 finished in 6.3 seconds.
+Chain 4 finished in 6.1 seconds.
 
 All 4 chains finished successfully.
-Mean chain execution time: 6.4 seconds.
-Total execution time: 25.9 seconds.</code></pre>
+Mean chain execution time: 6.3 seconds.
+Total execution time: 25.4 seconds.</code></pre>
 </div>
 <div class="sourceCode cell-code" id="cb74"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb74-1"><a href="#cb74-1" aria-hidden="true" tabindex="-1"></a>fit_mac_weibull_ind</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
 <div class="cell-output cell-output-stdout">