Summarizing and Comparing Models

Using bbr.bayes functions to summarize modeling history

bbr
model history
model comparison

Introduction

Summarizing and comparing multiple models in a workflow is often a relatively slow manual process. bbr has been extended to include functions that automate much of that process for Stan models. The resulting model summaries provide a condensed overview of the key decisions made during the model development process. The modeling summary functions take advantage of the annotation features of bbr (e.g., the description, notes, and tags discussed in the “Initial model” page) to summarize the important features of key models.

The page demonstrates how to:

  • Create and modify model run logs
  • Check and compare model tags
  • Calculate expected log predictive density (ELPD) for model comparison based on approximate leave one observation out and leave one individual (or group) out cross validation (LOO-CV and LOGO-CV)
  • Add LOO-CV ELPD and LOGO-CV ELPD to the summary tables constructed with bbr

Tools used

MetrumRG packages

bbr Manage, track, and report modeling activities, through simple R objects, with a user-friendly interface between R and NONMEM®.

bbr.bayes Extension of the bbr package to support Bayesian modeling with Stan or NONMEM®.

CRAN packages

dplyr A grammar of data manipulation.

Set up

Load required packages and set file paths to your model and figure directories.

library(bbr)
library(bbr.bayes)
library(here)
library(tidyverse)
library(yspec)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(tidybayes)
library(loo)
library(glue)
library(kableExtra)

set_cmdstan_path(path = here("Torsten", "v0.91.0", "cmdstan"))

MODEL_DIR <- here("model", "stan")
FIGURE_DIR <- here("deliv", "figure")

Using run_log and stan_summary_log

We can use run_log() to create a log of all, or a subset of, the models we’ve constructed.

model_log <- run_log(MODEL_DIR,
                     .include = paste0("ppkexpo", 1:5))

The return object is a tibble with information about each model including: (1) where the model resides on disk; (2) the model type (Stan vs stan_gq); (3) model description; (4) tags and notes added to the model object; (5) and a field indicating the model on which each model is based. For our worked example, a glimpse of model_log looks like this:

glimpse(model_log)
. Rows: 5
. Columns: 10
. $ absolute_model_path <chr> "/data/expo4-torsten-poppk/model/stan/ppkexpo1", "…
. $ run                 <chr> "ppkexpo1", "ppkexpo2", "ppkexpo3", "ppkexpo4", "p…
. $ yaml_md5            <chr> "1af6859e6361a6224be400c1dfcd4842", "9c800a3812b6a…
. $ model_type          <chr> "stan", "stan", "stan", "stan", "stan"
. $ description         <chr> "Base popPK model", "Base popPK model with non-cen…
. $ bbi_args            <list> <NULL>, <NULL>, <NULL>, <NULL>, <NULL>
. $ based_on            <list> <NULL>, "ppkexpo1", "ppkexpo2", "ppkexpo3", "ppkex…
. $ tags                <list> <"linear two compartment + first-order absorption…
. $ notes               <list> "Initial assessment of convergence based on R-hat…
. $ star                <lgl> FALSE, FALSE, FALSE, FALSE, FALSE

stan_add_summary() may be used to add information to the run log about the modeling results. This includes information about the Markov Chain Monte Carlo (MCMC) sampling and its performance and summary statistics for model parameters specified via the variables argument.

model_log <- model_log %>% 
  stan_add_summary(variables = c("lp__", "CLHat", "QHat", "V2Hat", "V3Hat",
                                "kaHat", "sigma", glue("omega[{1:5}]")),
                    summary_fns = list("median", "rhat"))

glimpse(model_log)
. Rows: 5
. Columns: 42
. $ absolute_model_path <chr> "/data/expo4-torsten-poppk/model/stan/ppkexpo1", "…
. $ run                 <chr> "ppkexpo1", "ppkexpo2", "ppkexpo3", "ppkexpo4", "p…
. $ yaml_md5            <chr> "1af6859e6361a6224be400c1dfcd4842", "9c800a3812b6a…
. $ model_type          <chr> "stan", "stan", "stan", "stan", "stan"
. $ description         <chr> "Base popPK model", "Base popPK model with non-cen…
. $ bbi_args            <list> <NULL>, <NULL>, <NULL>, <NULL>, <NULL>
. $ based_on            <list> <NULL>, "ppkexpo1", "ppkexpo2", "ppkexpo3", "ppkex…
. $ tags                <list> <"linear two compartment + first-order absorption…
. $ notes               <list> "Initial assessment of convergence based on R-hat…
. $ star                <lgl> FALSE, FALSE, FALSE, FALSE, FALSE
. $ method              <chr> "sample", "sample", "sample", "sample", "sample"
. $ stan_version        <chr> "2.33.0", "2.33.0", "2.33.0", "2.33.0", "2.33.0"
. $ nchains             <int> 4, 4, 4, 4, 4
. $ iter_warmup         <int> 500, 500, 500, 500, 500
. $ iter_sampling       <int> 500, 500, 500, 500, 500
. $ num_divergent       <int> 0, 0, 0, 0, 0
. $ num_max_treedepth   <int> 0, 0, 0, 0, 0
. $ bad_ebfmi           <lgl> TRUE, FALSE, FALSE, FALSE, FALSE
. $ lp___median         <dbl> 1217.71000, -31.81210, -31.61685, -30.02115, -30.…
. $ CLHat_median        <dbl> 3.108490, 3.141785, 3.133405, 3.278165, 3.273925
. $ QHat_median         <dbl> 3.798445, 3.741705, 3.769575, 3.766750, 3.754470
. $ V2Hat_median        <dbl> 60.84365, 61.63940, 61.77600, 61.74010, 61.65075
. $ V3Hat_median        <dbl> 67.71530, 67.49960, 67.76195, 67.80675, 67.69595
. $ kaHat_median        <dbl> 1.600785, 1.638600, 1.630025, 1.622190, 1.621595
. $ sigma_median        <dbl> 0.213561, 0.213678, 0.213322, 0.213414, 0.213452
. $ `omega[1]_median`   <dbl> 0.4351580, 0.4374375, 0.4065635, 0.3437790, 0.341…
. $ `omega[2]_median`   <dbl> 0.12347200, 0.11810250, 0.07907415, 0.05384745, 0.…
. $ `omega[3]_median`   <dbl> 0.3762690, 0.3767190, 0.2900915, 0.2888355, 0.2885…
. $ `omega[4]_median`   <dbl> 0.18236050, 0.18020650, 0.03654565, 0.03456670, 0.…
. $ `omega[5]_median`   <dbl> 0.5433645, 0.5364520, 0.5264730, 0.5178915, 0.5175…
. $ lp___rhat           <dbl> 1.3181695, 1.0081721, 0.9997599, 1.0053176, 1.0063…
. $ CLHat_rhat          <dbl> 1.001544, 1.082980, 1.071287, 1.055304, 1.024410
. $ QHat_rhat           <dbl> 1.015772, 1.004280, 1.010135, 1.001177, 1.003198
. $ V2Hat_rhat          <dbl> 1.002430, 1.026908, 1.039901, 1.014522, 1.008564
. $ V3Hat_rhat          <dbl> 1.0046124, 1.0042703, 1.0016776, 0.9993586, 1.0015…
. $ kaHat_rhat          <dbl> 1.003084, 1.008781, 1.013092, 1.005728, 1.000681
. $ sigma_rhat          <dbl> 1.002482, 1.001968, 1.002646, 1.001186, 1.001770
. $ `omega[1]_rhat`     <dbl> 1.003065, 1.046383, 1.011892, 1.008224, 1.004589
. $ `omega[2]_rhat`     <dbl> 1.165550, 1.012651, 1.008205, 1.006052, 1.003228
. $ `omega[3]_rhat`     <dbl> 1.010701, 1.007661, 1.013924, 1.006572, 1.006364
. $ `omega[4]_rhat`     <dbl> 1.0527630, 1.0013262, 1.0099005, 1.0041721, 0.9999…
. $ `omega[5]_rhat`     <dbl> 1.014980, 1.005679, 1.006394, 1.004728, 1.001840

From that tibble, it’s relatively easy to construct a table summarizing your model development. For example, the following code generates a table summarizing the characteristics of our sequence of models along with the posterior median log probability (lp__) for each model which is a crude measure of model fit.

model_log %>% 
  add_tags_diff() %>%
  select(run, based_on,  description, tags_added, tags_removed, lp___median) %>%
  mutate(lp___median = pmtables::sig(lp___median, 4)) %>%
  knitr::kable() %>%
  kable_styling()
run based_on description tags_added tags_removed lp___median
ppkexpo1 NULL Base popPK model linear two compartment + first-order absorption, centered parameterization , log(CL) ~ normal IIV , log(Q) ~ normal IIV , log(V2) ~ normal IIV , log(V3) ~ normal IIV , log(ka) ~ normal IIV , lognormal RUV 1218
ppkexpo2 ppkexpo1 Base popPK model with non-centered parameterization non-centered parameterization centered parameterization -31.81
ppkexpo3 ppkexpo2 PopPK model: ppkexpo2 + allometric scaling fixed allometric scaling -31.62
ppkexpo4 ppkexpo3 PopPK model: ppkexpo3 + allometric scaling & effects of EGFR, age, and albumin in CL EGFR on CL , age on CL , albumin on CL -30.02
ppkexpo5 ppkexpo4 PopPK model: ppkexpo4 without simulations for PPCs -30.44

Use stan_summary_log() to summarize and compare model results. The tibble it produces includes information about the MCMC sampling and its performance. It contains summary statistics for model parameters specified via the variables argument. Additionally, it contains the CmdStanMCMC object for each model.

sum_log <- stan_summary_log(MODEL_DIR,
                            .include = paste0("ppkexpo", 1:4))
glimpse(sum_log)
. Rows: 4
. Columns: 18
. $ absolute_model_path <chr> "/data/expo4-torsten-poppk/model/stan/ppkexpo1", "…
. $ run                 <chr> "ppkexpo1", "ppkexpo2", "ppkexpo3", "ppkexpo4"
. $ method              <chr> "sample", "sample", "sample", "sample"
. $ stan_version        <chr> "2.33.0", "2.33.0", "2.33.0", "2.33.0"
. $ nchains             <int> 4, 4, 4, 4
. $ iter_warmup         <int> 500, 500, 500, 500
. $ iter_sampling       <int> 500, 500, 500, 500
. $ num_divergent       <int> 0, 0, 0, 0
. $ num_max_treedepth   <int> 0, 0, 0, 0
. $ bad_ebfmi           <lgl> TRUE, FALSE, FALSE, FALSE
. $ lp___mean           <dbl> 1222.77933, -32.94283, -32.17479, -29.95754
. $ lp___median         <dbl> 1217.71000, -31.81210, -31.61685, -30.02115
. $ lp___q5             <dbl> 1066.71000, -82.62959, -77.54488, -74.23297
. $ lp___q95            <dbl> 1395.93950, 12.06736, 11.77985, 13.26280
. $ lp___rhat           <dbl> 1.3181695, 1.0081721, 0.9997599, 1.0053176
. $ lp___ess_bulk       <dbl> 10.25029, 294.79398, 407.64672, 397.98073
. $ lp___ess_tail       <dbl> 38.46279, 674.06085, 984.47737, 794.08719
. $ fit                 <list> <<CmdStanMCMC>
.   Inherits from: <CmdStanFit>
.   Publi…

The functions provided by the CmdStanMCMC objects in the summary log may be used to generate additional model results. For example, we demonstrate the use of the $loo() function to add information relevant to model comparison for the summary log. Particularly, we calculate the ELPD (commonly used metric for Bayesian (Bayes) model comparison) based on an approximation of LOO-CV.

sum_log <- sum_log %>%
  mutate(LOO_ELPD = map_vec(fit, \(x) x$loo(cores = 4, 
                                            save_psis = TRUE)$estimates[1, 1]))

sum_log %>%
  select(run, nchains, iter_warmup, iter_sampling, num_divergent, 
         num_max_treedepth, bad_ebfmi, LOO_ELPD) %>%
  knitr::kable() %>%
  kable_styling()
run nchains iter_warmup iter_sampling num_divergent num_max_treedepth bad_ebfmi LOO_ELPD
ppkexpo1 4 500 500 0 0 TRUE -18242.47
ppkexpo2 4 500 500 0 0 FALSE -18243.25
ppkexpo3 4 500 500 0 0 FALSE -18219.46
ppkexpo4 4 500 500 0 0 FALSE -18211.03

Model comparison using ELPD

LOO-CV approximates out of sample prediction error for a new observation in the same individuals. For mixed effects models like the population pharmacokinetic (popPK) models in this expo, a better metric for model comparison is based on LOGO-CV. In our case, the “group” of interest is an individual subject or patient. LOGO-CV approximates out of sample prediction error for new observations in new individuals.

The calculation of LOO-ELPD shown above uses the likelihoods contributed by each observation. For models ppkexpo1, ppkexpo2, ppkexpo3, and ppkexpo4, those likelihoods are contained in the log_lik vector. To calculate LOGO-ELPD, we need to derive the log-likelihood contribution for each independent unit. For our analysis, this corresponds to the likelihood at the subject level.

For non-linear mixed effect models, such as the one we’re using here, the log-likelihood is not available in a closed-form solution. So, we derive it using a Monte Carlo approximation. Specifically, we want to calculate the likelihood for the \(i^{th}\) individual:

\[ \log \ell_i(\boldsymbol{\theta},\boldsymbol{\Omega} ~|~ \mathbf{y}_i) = \log \left( \int f(\mathbf{y}_i ~|~ \boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\eta}) ~f(\boldsymbol{\eta} ~|~ \boldsymbol{\Omega}) ~ d\boldsymbol{\eta} \right) \] The Monte Carlo approximation to the integral is derived as:

\[\begin{align*} \ell_i(\boldsymbol{\theta},\boldsymbol{\Omega} ~|~ \mathbf{y}_i) &= \int \prod_{j=1}^{n_i} f(y_{ij} ~|~ \boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\eta}) ~f(\boldsymbol{\eta} ~|~ \boldsymbol{\Omega}) ~ d\boldsymbol{\eta} \\ & \approx \frac{1}{M} \sum_{m=1}^M \prod_{j=1}^{n_i} f(y_{ij} ~|~ \boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\eta}^{(m)}) ~f(\boldsymbol{\eta}^{(m)} ~|~ \boldsymbol{\Omega}) \\ & = \frac{1}{M} \sum_{m=1}^M \exp \left( \sum_{j=1}^{n_i} \log f(y_{ij} ~|~ \boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\eta}^{(m)}) + \log f(\boldsymbol{\eta}^{(m)} ~|~ \boldsymbol{\Omega}) \right) \end{align*}\]

Thus, the approximation to \(\log \ell_i(\boldsymbol{\theta},\boldsymbol{\Omega} ~|~ \mathbf{y}_i)\) is given by

\[\begin{align*} \log \ell_i(\boldsymbol{\theta},\boldsymbol{\Omega} ~|~ \mathbf{y}_i) & \approx \log(\sum_{m=1}^M \exp(A_i)) - \log M \\ & = \text{log\_sum\_exp}(A_i) - \log M \end{align*}\]

where \(M\) is the number of Monte Carlo samples, \(\eta^{(m)}\) is a draw from the distribution \(f(\boldsymbol{\eta} ~|~ \boldsymbol{\Omega})\), log_sum_exp corresponding to the Stan function of the same name, and

\[ A_i = \sum_{j=1}^{n_i} \log f(y_{ij} ~|~ \boldsymbol{\theta}, \boldsymbol{\Omega}, \boldsymbol{\eta}^{(m)}) + \log f(\boldsymbol{\eta}^{(m)} ~|~ \boldsymbol{\Omega}) \]

To calculate the LOGO-ELPD for our models, we take advantage of the ability in cmdstan to run a stand-alone generated quantities model like we did in the section on “Generated Quantities Models.” We demonstrate the process in detail with the ppkexpo2 model. Let’s start by re-reading the parent model and copying it to a generated quantities model called ppkexpo2_gq1.

ppkexpo2 <- read_model(here(MODEL_DIR, "ppkexpo2"))
ppkexpo2_gq1 <- copy_model_as_stan_gq(ppkexpo2, "ppkexpo2_gq1", .inherit_tags = TRUE, 
                      .overwrite = TRUE) %>%
  add_description("PopPK model: ppkexpo2 model to generate approximate individual log-likelihoods")

Now, manually edit ppkexpo2_gq1.stan by deleting the contents of the model{} block and adding the simulation code in the generated quantities{} block. Also edit ppkexpo2_gq1-standata.R to add idObs and idBlq to the data set. For this demo, we copy previously created files from the demo folder.

ppkexpo2_gq1 <- ppkexpo2_gq1 %>%
  add_stanmod_file(here(MODEL_DIR, "demo", "ppkexpo2_gq1.stan")) %>%
  add_standata_file(here(MODEL_DIR, "demo", "ppkexpo2_gq1-standata.R"))

Here is the resulting Stan model.

. /*
. Base PPK model
.   - Linear 2 compartment
.   - Non-centered parameterization
.   - lognormal residual variation
.   
. Based on NONMEM PopPK FOCE example project
. */
. 
. data{
.   int<lower = 1> nId; // number of individuals
.   int<lower = 1> nt; // number of events (rows in data set)
.   int<lower = 1> nObs; // number of PK observations
.   array[nObs] int<lower = 1> iObs; // event indices for PK observations
.   int<lower = 1> nBlq; // number of BLQ observations
.   array[nBlq] int<lower = 1> iBlq; // event indices for BLQ observations
.   array[nt] real<lower = 0> amt, time;
.   array[nt] int<lower = 1> cmt;
.   array[nt] int<lower = 0> evid;
.   array[nId] int<lower = 1> start, end; // indices of first & last observations
.   vector<lower = 0>[nObs] cObs;
.   array[nObs] int<lower = 1> idObs; // subject index associated with each cObs
.   array[nBlq] int<lower = 1> idBlq; // subject index associated with each BLQ
.   real<lower = 0> LOQ;
. }
. 
. transformed data{
.   int<lower = 1> nRandom = 5, nCmt = 3;
.   array[nt] real<lower = 0> rate = rep_array(0.0, nt), 
.                             ii = rep_array(0.0, nt);
.   array[nt] int<lower = 0> addl = rep_array(0, nt), 
.                            ss = rep_array(0, nt);
.   array[nId, nCmt] real F = rep_array(1.0, nId, nCmt), 
.                         tLag = rep_array(0.0, nId, nCmt);
.   }
. 
. parameters{
.   real<lower = 0> CLHat, QHat, V2Hat, V3Hat;
. // To constrain kaHat > lambda1 uncomment the following
.   real<lower = (CLHat / V2Hat + QHat / V2Hat + QHat / V3Hat +
.       sqrt((CLHat / V2Hat + QHat / V2Hat + QHat / V3Hat)^2 -
.            4 * CLHat / V2Hat * QHat / V3Hat)) / 2> kaHat; // ka > lambda_1
.   
.   real<lower = 0> sigma;
.   vector<lower = 0>[nRandom] omega;
. //  corr_matrix[nRandom] rho;
. //  array[nId] vector[nRandom] theta;
.   cholesky_factor_corr[nRandom] L_corr;
.   matrix[nRandom, nId] z;
. }
. 
. transformed parameters{
.   vector[nRandom] thetaHat = log([CLHat, QHat, V2Hat, V3Hat, kaHat]');
.   // Individual parameters
.   matrix[nId, nRandom] theta = (rep_matrix(thetaHat, nId) + diag_pre_multiply(omega, L_corr * z))';
.   vector<lower = 0>[nId] CL = exp(theta[,1]),
.                          Q = exp(theta[,2]),
.                          V2 = exp(theta[,3]),
.                          V3 = exp(theta[,4]),
.                          ka = exp(theta[,5]);
.   corr_matrix[nRandom] rho =  L_corr * L_corr';
.   cov_matrix[nRandom] Omega = quad_form_diag(rho, omega);
. 
.   row_vector[nt] cHat; // predicted concentration
.   matrix[nCmt, nt] x; // mass in all compartments
.   
.   for(j in 1:nId){
.     x[, start[j]:end[j]] = pmx_solve_twocpt(time[start[j]:end[j]],
.                                             amt[start[j]:end[j]],
.                                             rate[start[j]:end[j]],
.                                             ii[start[j]:end[j]],
.                                             evid[start[j]:end[j]],
.                                             cmt[start[j]:end[j]],
.                                             addl[start[j]:end[j]],
.                                             ss[start[j]:end[j]],
.                                             {CL[j], Q[j], V2[j], V3[j], ka[j]});
. 
.     cHat[start[j]:end[j]] = x[2, start[j]:end[j]] / V2[j];
.   }
. }
. 
. model{
. 
. }
. 
. generated quantities{
.   // approximate individual log-likelihood
.   vector[nId] log_lik_id;
.   {
.     int nSim = 1000;
.     // simulated individual log-likelihood
.     array[nId] vector[nSim] log_lik_id_sim;
.     // Individual parameters for hypothetical new individuals
.     array[nId] vector[nRandom] thetaNew;
.     array[nId] real CLNew, QNew, V2New, V3New, kaNew;
.     
.     row_vector[nt] cHatNew; // predicted concentration in new individuals
.     matrix[nCmt, nt] xNew; // mass in all compartments in new individuals
. 
.     for(m in 1:nSim){
.       thetaNew = multi_normal_rng(rep_array(thetaHat, nId), Omega);
.       CLNew = exp(thetaNew[,1]);
.       QNew = exp(thetaNew[,2]);
.       V2New = exp(thetaNew[,3]);
.       V3New = exp(thetaNew[,4]);
.       kaNew = exp(thetaNew[,5]);
.       for(j in 1:nId){
.         xNew[, start[j]:end[j]] = pmx_solve_twocpt(time[start[j]:end[j]],
.         amt[start[j]:end[j]],
.         rate[start[j]:end[j]],
.         ii[start[j]:end[j]],
.         evid[start[j]:end[j]],
.         cmt[start[j]:end[j]],
.         addl[start[j]:end[j]],
.         ss[start[j]:end[j]],
.         {CLNew[j], QNew[j], V2New[j], 
.         V3New[j], kaNew[j]});
.         
.         cHatNew[start[j]:end[j]] = xNew[2, start[j]:end[j]] / V2New[j];
.       }
.       
.       // Calculate log-likelihood values to use for LOO CV
.       
.       // Individual-level log-likelihood
.       for(j in 1:nId)
.         log_lik_id_sim[j, m] = multi_normal_lpdf(thetaNew[j] | thetaHat, Omega);
.       for(i in 1:nObs)
.         log_lik_id_sim[idObs[i], m] += lognormal_lpdf(cObs[i] | log(cHatNew[iObs[i]]), sigma);
.       for(i in 1:nBlq)
.         log_lik_id_sim[idBlq[i], m] += lognormal_lcdf(LOQ | log(cHatNew[iBlq[i]]), sigma);
.     }
.     for(j in 1:nId){
.       log_lik_id[j] = log_sum_exp(log_lik_id_sim[j,]) - log(nSim);
.     }
.   }
. }

Submit the model

## Set cmdstanr arguments
ppkexpo2_gq1 <- ppkexpo2_gq1 %>%
  set_stanargs(list(parallel_chains = 4,
                    seed = 1234),
               .clear = TRUE)

## Check that the necessary files are present
check_stan_model(ppkexpo2_gq1)

## Simulate using Stan
ppkexpo2_gq1_sim <- ppkexpo2_gq1 %>% submit_model(.overwrite = TRUE)

ppkexpo2_loo_id <- loo(ppkexpo2_gq1_sim$draws(variables = "log_lik_id"), cores = 4, save_psis = TRUE)
ppkexpo2_loo_id

Now, let’s do the same for ppkexpo3_gq1 and ppkexpo4_gq1.

ppkexpo3 <- read_model(here(MODEL_DIR, "ppkexpo3"))
ppkexpo3_gq1 <- copy_model_as_stan_gq(ppkexpo3, "ppkexpo3_gq1", .inherit_tags = TRUE, 
                      .overwrite = TRUE) %>%
  add_description("PopPK model: ppkexpo3 model to generate approximate individual log-likelihoods")
ppkexpo3_gq1 <- ppkexpo3_gq1 %>%
  add_stanmod_file(here(MODEL_DIR, "demo", "ppkexpo3_gq1.stan")) %>%
  add_standata_file(here(MODEL_DIR, "demo", "ppkexpo3_gq1-standata.R"))

ppkexpo4 <- read_model(here(MODEL_DIR, "ppkexpo4"))
ppkexpo4_gq1 <- copy_model_as_stan_gq(ppkexpo4, "ppkexpo4_gq1", .inherit_tags = TRUE, 
                      .overwrite = TRUE) %>%
  add_description("PopPK model: ppkexpo4 model to generate approximate individual log-likelihoods")
ppkexpo4_gq1 <- ppkexpo4_gq1 %>%
  add_stanmod_file(here(MODEL_DIR, "demo", "ppkexpo4_gq1.stan")) %>%
  add_standata_file(here(MODEL_DIR, "demo", "ppkexpo4_gq1-standata.R"))

ppkexpo3_gq1 <- ppkexpo3_gq1 %>%
  set_stanargs(list(parallel_chains = 4,
                    seed = 1234),
               .clear = TRUE)

## Check that the necessary files are present
check_stan_model(ppkexpo3_gq1)

## Simulate using Stan
ppkexpo3_gq1_sim <- ppkexpo3_gq1 %>% submit_model(.overwrite = TRUE)

ppkexpo4_gq1 <- ppkexpo4_gq1 %>%
  set_stanargs(list(parallel_chains = 4,
                    seed = 1234),
               .clear = TRUE)

## Check that the necessary files are present
check_stan_model(ppkexpo4_gq1)

## Simulate using Stan
ppkexpo4_gq1_sim <- ppkexpo4_gq1 %>% submit_model(.overwrite = TRUE)

If ppkexpo2, ppkexpo3, ppkexpo4, ppkexpo2_gq1, ppkexpo3_gq1, and ppkexpo4_gq1 have already been run, you can read them using the following statements without having to re-simulate.

ppkexpo2 <- read_model(here(MODEL_DIR, "ppkexpo2"))
ppkexpo2_fit <- read_fit_model(ppkexpo2)

ppkexpo3 <- read_model(here(MODEL_DIR, "ppkexpo3"))
ppkexpo3_fit <- read_fit_model(ppkexpo3)

ppkexpo4 <- read_model(here(MODEL_DIR, "ppkexpo4"))
ppkexpo4_fit <- read_fit_model(ppkexpo4)

ppkexpo2_gq1 <- read_model(here(MODEL_DIR, "ppkexpo2_gq1"))
ppkexpo2_gq1_sim <- read_fit_model(ppkexpo2_gq1)

ppkexpo3_gq1 <- read_model(here(MODEL_DIR, "ppkexpo3_gq1"))
ppkexpo3_gq1_sim <- read_fit_model(ppkexpo3_gq1)

ppkexpo4_gq1 <- read_model(here(MODEL_DIR, "ppkexpo4_gq1"))
ppkexpo4_gq1_sim <- read_fit_model(ppkexpo4_gq1)

Calculate LOO-ELPD and LOGO-ELPD

Now, let’s calculate LOO-ELPD and LOGO-ELPD for each model. As before, we calculate LOO-ELPD using the loo function of the fit objects for each model. For LOGO-ELPD, we need to call the loo function of the loo R package.

sum_log <- stan_summary_log(MODEL_DIR,
                            .include = paste0("ppkexpo", 2:4))


## LOO-CV 
ppkexpo2_loo <- ppkexpo2_fit$loo(cores = 4, save_psis = TRUE)
ppkexpo3_loo <- ppkexpo3_fit$loo(cores = 4, save_psis = TRUE)
ppkexpo4_loo <- ppkexpo4_fit$loo(cores = 4, save_psis = TRUE)

## LOGO-CV 
ppkexpo2_loo_id <- loo(ppkexpo2_gq1_sim$draws(variable = "log_lik_id"), cores = 4, save_psis = TRUE)
ppkexpo3_loo_id <- loo(ppkexpo3_gq1_sim$draws(variable = "log_lik_id"), cores = 4, save_psis = TRUE)
ppkexpo4_loo_id <- loo(ppkexpo4_gq1_sim$draws(variable = "log_lik_id"), cores = 4, save_psis = TRUE)

Add LOO-ELPD and LOGO-ELPD to the summary log

Now, we can add both LOO-ELPD and LOGO-ELPD to the summary log and print it as a table.

loo_id <- data.frame(run = c("ppkexpo2", "ppkexpo3", "ppkexpo4"),
                     elpd_loo = c(ppkexpo2_loo$estimates[1,1],
                                   ppkexpo3_loo$estimates[1,1],
                                   ppkexpo4_loo$estimates[1,1]),
                     elpd_logo = c(ppkexpo2_loo_id$estimates[1,1],
                                   ppkexpo3_loo_id$estimates[1,1],
                                   ppkexpo4_loo_id$estimates[1,1]))
sum_log <- sum_log %>%
  left_join(loo_id)

sum_log %>%
  select(run, iter_warmup, iter_sampling, num_divergent, bad_ebfmi, elpd_loo, elpd_logo) %>%
  mutate(elpd_loo = pmtables::sig(elpd_loo, 5),
         elpd_logo = pmtables::sig(elpd_logo, 5)) %>%
  knitr::kable() %>%
  kable_styling()
run iter_warmup iter_sampling num_divergent bad_ebfmi elpd_loo elpd_logo
ppkexpo2 500 500 0 FALSE -18243 -18474
ppkexpo3 500 500 0 FALSE -18219 -18075
ppkexpo4 500 500 0 FALSE -18211 -17938

Compare ppkexpo2, ppkexpo3, and ppkexpo4

The loo_compare function calculates various quantities relevant to comparing the models. These include the difference between the ELPD value for each model relative to the model with the highest value and the standard error of that difference. Let’s apply loo_compare to both LOO-ELPD and LOGO-ELPD.

Model comparison with respect to LOO-CV

loo_compare(ppkexpo2_loo, ppkexpo3_loo, ppkexpo4_loo) %>%
  knitr::kable(digits = 2) %>%
  add_header_above(c("model1 = ppkexpo2, model2 = ppkexpo3, model3 = ppkexpo4" = 9)) %>%
  add_header_above(c("LOO-CV" = 9), bold = TRUE) %>%
  kable_styling()
LOO-CV
model1 = ppkexpo2, model2 = ppkexpo3, model3 = ppkexpo4
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
model3 0.00 0.00 -18211.03 99.83 336.16 12.32 36422.06 199.67
model2 -8.43 4.99 -18219.46 99.82 347.49 12.55 36438.92 199.65
model1 -32.22 8.10 -18243.25 99.16 379.75 13.04 36486.50 198.32

Model comparison with respect to LOGO-CV

loo_compare(ppkexpo2_loo_id, ppkexpo3_loo_id, ppkexpo4_loo_id) %>%
  knitr::kable(digits = 2) %>%
  add_header_above(c("model1 = ppkexpo2, model2 = ppkexpo3, model3 = ppkexpo4" = 9)) %>%
  add_header_above(c("LOGO-CV" = 9), bold = TRUE) %>%
  kable_styling()
LOGO-CV
model1 = ppkexpo2, model2 = ppkexpo3, model3 = ppkexpo4
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
model3 0.00 0.00 -17938.07 538.69 705.38 6.84 35876.14 1077.39
model2 -137.21 23.07 -18075.28 542.79 504.27 9.04 36150.56 1085.58
model1 -535.80 28.95 -18473.86 543.99 263.29 13.28 36947.73 1087.98

The LOGO-ELPD indicates that ppkexpo4 results in a greater LOGO-ELPD than the other models, and the difference is substantially greater than its standard error. The LOO-ELPD results are less definitive, but that’s not too surprising given that the models differ primarily with respect to model elements related to interindividual variation. Such differences are expected to have greater effect on out of sample error for new individuals than that for just new observations in the same individuals.

Other resources

The following script from the GitHub repository is discussed on this page. If you’re interested running this code, visit the About the GitHub Repo page first.

Modeling Summary script: modeling-summary.Rmd