In mhismail/mrgsolvetk: Simulation Tool Kit for mrgsolve

knitr::opts_chunk$set(comment='.', message=FALSE, 
                      fig.path="inst/maintenance/img/README-")

mrgsolvetk

A toolkit to be used with mrgsolve

Installation

library(devtools)
install_github("mrgsolve/mrgsolvetk", ref = "mrgoptim")

Examples

• Sensitivity Analyses
• Parameter Optimization

library(ggplot2)
library(dplyr)
library(mrgsolve)
library(mrgsolvetk)

theme_set(theme_bw())

mod <- mread_cache("pk1cmt",modlib())

mod <- ev(mod, amt=100) %>% Req(CP) %>% update(end = 48, delta = 0.25)

param(mod)

Sensitivity analyses

sens_unif

• Draw parameters from uniform distribution based on current parameter values
• lower and upper scale the parameter value to provide a and b arguments to runif

out <- 
  mod %>% 
  select(CL,VC,KA1) %>%
  sens_unif(.n=10, lower=0.2, upper=3)

out

sens_plot(out, CP)

We can also make a univariate version of this

mod %>% 
  select(CL,VC,KA1) %>%
  sens_unif(.n=10, lower=0.2, upper=3, univariate = TRUE) %>%
  sens_plot(CP, split = TRUE)

sens_norm

• Draw parameters from (log) normal distribution based on current parameter values and %CV

mod %>% 
  select(CL,VC) %>%
  sens_norm(.n=10, cv=30) %>%
  sens_plot(CP)

sens_seq

• Give a sequence for one or more parameters

mod %>% sens_seq(CL = seq(2,12,2), VC = seq(30,100,10)) %>% sens_plot(CP)

sens_range

• Create sets of parameters equally-spaced between two bounds

mod %>%
  select(CL,VC) %>%
  sens_range(.n = 5, .factor = 4) %>%
  sens_plot(CP, split = TRUE)

or

mod %>%
  sens_range(CL = c(0.5, 1.5), VC = c(10,40), .n = 5) %>%
  sens_plot(CP)

sens_grid

• Like sens_seq but performs all combinations

mod %>%  sens_grid(CL = seq(1,10,1), VC = seq(20,40,5)) %>% sens_plot(CP)

sens_covset

• Use dmutate to generate random variates for each parameter

cov1 <- dmutate::covset(CL ~ runif(1,3.5), VC[0,] ~ rnorm(50,25))

cov1

out <- mod %>% sens_covset(cov1) 

out

distinct(out,ID,CL,VC)

Maximum Likelihood Parameter Optimization

mrgoptim

This example shows a simultaneous fit of PK and PD data from five dose levels.

Data structure

The data to be fit is an mrgsolve dataset. Required columns for fitting are:

• ID
• time
• evid
• cmt
• amt
• dv

data <- read.csv("inst/maintenance/data/optim-example.csv")

head(data)

Plot the data to get an idea of the profiles to be fit. cmt 1 is plasma concentration data and cmt 2 is PD data

ggplot(data, aes(x = time, y = dv, color = as.factor(ID))) +
  geom_point() +
  geom_line() +
  facet_wrap("cmt") +
  guides(color = FALSE)

The following model will be fit to these data:

• PK: 2 compartment model
• PD: Inhibitory Emax model with baseline
• Proportional error models for both PK and PD

code<-"
$PROB 2 cmt PK Model, Emax PD model

$PARAM
CL=10
VC = 20
VP = 20
Q=20
Emax = 60
BL = 50
EC50 = 10
gamma =1
sigma1 = 0.1
sigma2 = 0.1

$CMT X1 X2 

$ODE
dxdt_X1 = -(Q+CL)/VC*X1+Q/VP*X2;
dxdt_X2 = Q/VC*X1-Q/VP*X2;

$TABLE
capture PK = X1/VC;
capture varPK = (PK*sigma1)*(PK*sigma1);


capture PD = BL-(pow(PK,gamma)*Emax)/(pow(PK,gamma)+pow(EC50,gamma));
capture varPD = (PD*sigma2)*(PD*sigma2);

capture ipred = NAN;
capture var = NAN; 


if(self.cmt == 1) {
   ipred = PK;
   var = varPK;
}

if(self.cmt == 2) {
   ipred = PD;
   var = varPD;
}"

mod <- mcode("2cmtPK-Emax", code)

Here, the predicted plasma concentrations, response, and variances were captured in the PK, PD, varPK, and varPD outputs, respectively. Predictions and variances are consolidated to a single column each. If cmt == 1 the predicted output, ipred, will be PK and prediction variance, var, varPK. If cmt == 2 the predicted output will be PD and prediction variance varPD.

Let's check how the initial parameter values fit the data.

out <- mod %>%
  data_set(data) %>%
  carry.out(cmt, dv) %>%
  obsonly() %>%
  mrgsim() %>%
  as.data.frame()

ggplot(filter(out, cmt == 1), aes(x = time, y = ipred, color = as.factor(ID))) +
  geom_line() +
  geom_point(aes(y = dv)) +
  guides(color = FALSE)

ggplot(filter(out, cmt == 2), aes(x = time, y = ipred, color = as.factor(ID))) +
  geom_line() +
  geom_point(aes(y = dv)) +
  guides(color = FALSE)

Not terrible, should be good enough for initial estimates.

Now let's use mrgoptim to optimize the parameters and return parameter values and precision. Use the output, and var arguments to specify which columns in the model code correspond to the predicted values and variances. Specify which system parameters to optimize with the prms argument and variance parameters with the v_prms arguments.

fit <- mod %>%
  data_set(data) %>%
  mrgoptim(output = "ipred",
           var = "var",
           prms = c("CL",
                    "VC",
                    "VP",
                    "Q",
                    "Emax",
                    "BL",
                    "EC50",
                    "gamma"),
           v_prms = c("sigma1", "sigma2"),
           method = "newuoa")

The function returns a list with some information about the optimization, the final objective function value (-LL), final parameter estimates, covariance and correlation matrices, CV percent, and output dataset.

print(fit)

Lets check how the optimized parameters fit the data.

out_fit <- mod %>%
  param(fit$par) %>%
  carry.out(cmt, dv) %>%
  data_set(data) %>%
  obsonly() %>%
  mrgsim() %>%
  as.data.frame()


ggplot(filter(out_fit, cmt == 1), aes(x = time, y = ipred, color = as.factor(ID))) +
  geom_line() +
  geom_point(aes(y = dv)) +
  guides(color = FALSE)

ggplot(filter(out_fit, cmt == 2), aes(x = time, y = ipred, color = as.factor(ID))) +
  geom_line() +
  geom_point(aes(y = dv)) +
  guides(color = FALSE)

mhismail/mrgsolvetk documentation built on May 7, 2023, 1:52 p.m.