In mrgsolve/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")

To install branch with mrgoptim function, until merged with mrgsolve/master:

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

Examples

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=F)

The following model will be fit to these data:

• PK: 2 compartment model
• PD: 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);"

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

Here, the plasma concentrations, response, and variances were captured in the PK, PD, varPK, and varPD outputs, respectively.

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

out <- mod%>%
  data_set(data)%>%
  obsonly()%>%
  mrgsim()%>%
  as.data.frame()

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

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

Not terrible, should be good enough for initial estimates.

Now let's use mrgoptim to optimize the parameters and return parameter values and precision. The input, output,and var arguments map the observed values to the output compartments and variances of the outputs. Since compartment 1 in the input dataset corresponds to plasma concentration data, which is the PK column in the output, they both need to be the first value in their respective vectors. Similarly, since varPK corresponds to the variance of the PK data, it is the first value in the var argument vector.

fit<- mod%>%
  data_set(data)%>%
  mrgoptim(input=c(1,2),
           output=c("PK","PD"),
           var= c("varPK","varPD"),
           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 (-2LL), 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)%>%
  data_set(data)%>%
  obsonly()%>%
  mrgsim()%>%
  as.data.frame()


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

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

mrgsolve/mrgsolvetk documentation built on May 11, 2019, 4:19 p.m.