Outline

1. Markov Chain Monte Carlo

• Principle & Workflow

• Example: linear model

2. MCMC()

• Working with GNU MCSim & R

• Example: Ethylbenzene PBPK model

3. Population modeling

• Single and multi level

• Example: Tetrachloroethylene PBPK model

General workflow

1 Model constructing or translating

2 Verify modeling result

• Compare with published result
• Mass balance

3 Uncertainty and sensitivity analysis

• Monte Carlo simulation
• Morris elementary effects screening
• Fourier amplitude sensitivity test

4 Model calibration

• Markov chain Monte Carlo
  • Diagnostics (Goodness-of-fit, convergence)

MC & MCMC

Monte Carlo

- Generate the probability distribution based on the given condition

Markov Chain Monte Carlo

- Iterative update the probability distribution, if the new proposed distribution is accepted.

Bayes' rule

$\theta$: Observed or unobserved parameter

$y$: Observed data

$p\left(\theta \right)$: Prior distribution of model parameter

$p\left(y|\theta \right)$: Sampling distribution of the experiment data given by a parameter vector

$p\left(y\right)$: Likelihood of data

$p\left(\theta |y\right)$: Posterior distribution

Monty Hall Problem

Suppose you’re on a game show, and you’re given the choice of three doors. Behind one door is a car, behind the others, goats. You pick a door, say number 1, and the host, who knows what’s behind the doors, opens another door, say number 3, which has a goat. He says to you, ”Do you want to pick door number 2?” Is it to your advantage to switch your choice of doors?

Marilyn vos Savant. 1990. Ask Marilyn. Parade Magazine: 16.

Markov Chain Monte Carlo

The algorithm was named for Nicholas Metropolis (physicist) and Wilfred Keith Hastings (statistician). The algorithm proceeds as follows.

Initialize

1. Pick an initial parameter sets ${\theta }_{t=0}=\{{\theta }_1,{\theta }_2,...{\theta }_n\}$

Iterate

1. Generate: randomly generate a candidate parameter state ${\theta }^{\prime }$ for the next sample by picking from the conditional distribution $J\left({\theta }^{\prime }|{\theta }_t\right)$
2. Compute: compute the acceptance probability $A({\theta }^{\prime },{\theta }_t)=min(1,\frac{P\left({\theta }^{\prime }|y\right)}{P\left({\theta }_t|y\right)}\frac{J\left({\theta }_t|{\theta }^{\prime }\right)}{J\left({\theta }^{\prime }|{\theta }_t\right)})$
3. Accept or Reject:
   1. generate a uniform random number $u\in [0,1]$
   2. if $u\le A(x^{\prime },x_t)$ accept the new state and set ${\theta }_{t+1}={\theta }^{\prime }$, otherwise reject the new state, and copy the old state forward ${\theta }_{t+1}={\theta }_t$

log-likelihood

The log-likelihood function was used to assess the goodness-of-fit of the model to the data (Woodruff and Bois, 1993, Hsieh et al., 2018)

$$LL=\sum _{i=1}^N-\frac{1}{2}\cdot \frac{{(y_i-{\hat{y}}_i)}^2}{S_{j\left[i\right]}^2}-\frac{1}{2}\ln \left(2\pi s_{j\left[i\right]}^2\right)$$

$N$: the total number of the data points

$y_i$: experimental observed

${\hat{y}}_i$: model predicted value

$j\left[i\right]$: data type (output variable) of data point $i$

$S_{j\left[i\right]}^2$: the variance for data type $j$

Calibration & evaluation

Prepare model and input files

• Need at least 4 chains in simulation

Check convergence & graph the output result

Evaluate the model fit

• Global evaluation (all output variables, $y_{i,t}$)
• Individual evaluation

Example: linear modeling

## linear.model.R ####
Outputs = {y}
# Model Parameters
A = 0; # Default value of intercept
B = 1; # Default value of slope
# Statistical parameter
SD_true = 0;
CalcOutputs { 
  y = A + B * t + NormalRandom(0,SD_true); 
}
End.

## ./mcsim.linear.model.R.exe linear_mcmc.in.R ####
MCMC ("MCMC.default.out","", # name of output file
     "",         # name of data file
     2000,0,     # iterations, print predictions flag
     1,2000,     # printing frequency, iters to print
     10101010);  # random seed (default)
Level {
  Distrib(A, Normal, 0, 2); # prior of intercept 
  Distrib(B, Normal, 1, 2); # prior of slope 
  Likelihood(y, Normal, Prediction(y), 0.5);
  Simulation {
    PrintStep (y, 0, 10, 1); 
    Data (y, 0.0, 0.15, 2.32, 4.33, 4.61, 6.68, 7.89, 7.13, 7.27, 9.4, 10.0);
  }
}
End.

The data

x <- seq(0, 10, 1)
y <- c(0.0, 0.15, 2.32, 4.33, 4.61, 6.68, 7.89, 7.13, 7.27, 9.4, 10.0)
plot(x, y)

MCMC sampling process

##   iter     A.1.     B.1.    LnPrior    LnData LnPosterior
## 1    0 -0.95484 -1.22704  -3.958099 -4567.137   -4571.095
## 2    1  1.83302 -1.22704  -4.264131 -3201.785   -3206.049
## 3    2  1.83302 -1.22704  -4.264131 -3201.785   -3206.049
## 4    3  4.78650 -1.22704  -6.707952 -2128.378   -2135.086
## 5    4  4.78650 -1.22704  -6.707952 -2128.378   -2135.086
## 6    5  7.38065 -1.22704 -10.653380 -1502.173   -1512.826

Trace plot

parms_name <- c("A.1.","B.1.")
mcmc_trace(sims, pars = parms_name, facet_args = list(ncol = 2))

Kernel density

j <- c(1002:2001) # Discard first half as burn-in
mcmc_dens_overlay(x = sims[j,,], pars = parms_name)

Pair plot

mcmc_pairs(sims[j,,], pars = parms_name, off_diag_fun = "hex")

Summary report

monitor(sims, probs = c(0.025, 0.975) , digit=4)

## Inference for the input samples (4 chains: each with iter=2001; warmup=1000):
## 
##                  mean se_mean       sd      2.5%     97.5% n_eff   Rhat
## iter        1500.0000 91.1594 288.9998 1025.0000 1975.0000    10 2.1056
## A.1.           0.4466  0.0300   0.3002   -0.1278    0.9969   100 1.0407
## B.1.           0.9980  0.0049   0.0512    0.8949    1.0919   109 1.0352
## LnPrior       -3.2607  0.0039   0.0377   -3.3492   -3.2247    94 1.0412
## LnData       -19.5507  0.0839   1.0927  -22.8720  -18.4954   170 1.0194
## LnPosterior  -22.8113  0.0868   1.1122  -26.0994  -21.7368   164 1.0203
## 
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Evaluation of model fit

Demo - linear model

1 Single chain testing

• MCMC simulation

2 Multi-chains simulation

• Check convergence

• Evaluation of model fit

MCMC()

# Input-file
MCMC();
# <Global assignments and specifications>
Level {
  Distrib();  
  Likelihood();
  # Up to 10 levels of hierarchy
  Simulation {
    # <Local assignments and specifications>
  }
  Simulation {
    # <Local assignments and specifications>
  }
  # Unlimited number of simulation specifications
} # end Level
End.

MCMC()

The statement, gives general directives for MCMC simulations with following syntax:

 MCMC("<OutputFilename>", "<RestartFilename>", "<DataFilename>",
          <nRuns>, <simTypeFlag>, <printFrequency>, <itersToPrint>,
          <RandomSeed>);

"<OutputFilename>" Output file name, the default is "MCMC.default.out"

"<RestartFilename>" Restart file name

"<DataFilename>" Data file name

<nRuns> an integer for the total sampling number (iteration)

<simTypeFlag> an integer (from 0 to 5) to define the simulation type

<printFrequency> an integer to set the interval of printed output

<itersToPrint> an integer to set the number of printed output from the final iteration

<RandomSeed> a numeric for pseudo-random number generator

Simulation types

<simTypeFlag> an integer (from 0 to 5) to define the simulation type

0, start/restart a new or unfinished MCMC simulations

1, use the last MCMC iteration to quickly check the model fit to the data

2, improve computational speed when convergence is approximately obtained

3, tempering MCMC with whole posterior

4, tempering MCMC with only likelihood

5, stochastic optimization

Check convergence

Manipulate (MCSim under R)

mcmc_array()

Visualize (bayesplot, corrplot)

mcmc_trace()

mcmc_dens_overlay()

mcmc_pairs()

corrplot()

Report (rstan)

monitor()

Parallel

Demo - Ethylbenzene PBPK Model

Model verification (tutorial 1)

- Compare the simulated result with previous study

Uncertainty analysis (tutorial 2)

- Set the probability distribution for model parameters

Morris elementary effects screening (tutorial 2)

- Find the influential parameters

MCMC calibration (tutorial 3)

- Estimate the "posterior"

Bayesian Population Model

Individuals level

$E$: Exposure

$t$: Time

$\phi$: measured parameters

$\theta$: unknown parameters

$y$: condition the data

Population level

$\mu$: Population means

${\Sigma }^2$: Population variances

${\sigma }^2$: Residual errors

Population modeling

# Comments of Input-file
<Global assignments and specifications>
   Level { # priors on population parameters
    Distrib ()
    Likelihood ()
     Level { # all subjects grouped
        Experiment { # 1 
          <Specifications and data in 1st simulation>
        } # end 1 
        Experiment { # 2
          <Specifications and data in 2nd simulation>
        } # end 2
      # ... put other experiments
    } # End grouped level  
  } # End population level
End.

# Comments of Input-file
<Global assignments and specifications>
   Level { # priors on population parameters
    Distrib ()
    Likelihood ()
     Level { # individuals
      Distrib ()
       Level { # subject A
        Experiment { # 1 
          <Specifications and data in 1st simulation>
        } # end 1 
        Experiment { # 2
          <Specifications and data in 2nd simulation>
        } # end 2
      } # End subject A
      # ... put other subjects
    } # End individuals level  
  } # End population level
End.

Demo - Tetrachloroethylene PBPK Model

1 Individuls level

• Population-experiments

2 Subject level

• Population-individuals-experiments

https://rpubs.com/Nanhung/SM_190523

Final thoughts

MCMC simulation is a computational intensive process for metamodel or big data, one of the important issue today is how to improve the efficiency in the modeling.

• If possible, remove or turn-off unnecessary state variable (e.g., mass balance, metabolites).

• Re-parameterization (Chiu et al., 2006)

• Global sensitivity analysis and parameter fixing (Hsieh et al., 2018)

• Parallel computing

• Grouped data

• Algorithm

• Software, hardware, ...