README.Rmd

---
output: github_document
title: ""
author: ""
date: ""
---
```{r,echo=FALSE}
knitr::opts_chunk$set(comment='.', message=FALSE, 
                      fig.path="inst/maintenance/img/README-")
```



# mrgsolvetk
A toolkit to be used with `mrgsolve`

# Installation
```{r, eval = FALSE}
library(devtools)
install_github("mrgsolve/mrgsolvetk", ref = "mrgoptim")
```


# Examples

- [Sensitivity Analyses](#sensitivity-analyses)
- [Parameter Optimization](#maximum-likelihood-parameter-optimization)
```{r warning=FALSE}
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`

```{r}
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

```{r}

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`

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

### `sens_seq`

- Give a sequence for one or more parameters

```{r}
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

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

or

```{r}

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

```{r}
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 

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

cov1
```

```{r}
out <- mod %>% sens_covset(cov1) 
```

```{r}
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	


```{r}
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

```{r warning=FALSE}
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

```{r}
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. 

```{r warning=FALSE}
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. 

```{r echo=T, results='hide'}
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.

```{r}
print(fit)

```
Lets check how the optimized parameters fit the data. 
```{r warning=FALSE}
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)


```