---
title: "Design evaluation and optimization in discrete space"
classoption: openany
output:
  rmarkdown::html_vignette:
    toc: true
bibliography: references.bib
biblio-style: apalike
link-citations: yes
linkcolor: blue
urlcolor: green
vignette: >
  %\VignetteIndexEntry{Design evaluation and optimization in discrete space}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

  <style type="text/css">
  
  body{
  font-size: 11pt;
  }
  
  h1.title {
  font-size: 40pt;
  }
  
  code.r{
  font-size: 11pt;
  }
  
  body {
  max-width: 1600px;  
  margin-left:10px;
  line-height: 20px;
  }
  
  </style>
  
```{r global_options, echo = FALSE, include = FALSE}
backup_options <- options()
options(width = 3000)
knitr::opts_chunk$set(collapse = TRUE,
comment = "#>",echo = FALSE, warning = FALSE, message = FALSE,
cache = FALSE, tidy = FALSE, size = "small")
library(PFIM)
```

# Overview

Data for this example result from a PKPD population study on a non-steroidal molecule [@FloresMurrieta1998]. Single doses of 1, 3.2, 10, 31.6, 56.2, or 100mg/kg per os were given respectively to 6 groups, each of which contained at least 6 rats (weighing 200g on average), following a parallel design. Blood sampling and drug response (DI score) evaluation were conducted at 0, 15, 30, and 45 min and at 1, 1.25, 1.5, 2, 3 and 4 hours after administration. Based on the evaluation results, we choose to optimize a design for 30 rats receiving a single dose of either 100mg/kg or 320 mg/kg, selecting only 3 from previously defined time points for blood sampling and response evaluation by using Fedorov-Wynn method and Multiplicative Algorithm. Reports of the design evaluation and optimization are available at <https://github.com/packagePFIM>

# Design evaluation 

### Define the model equations of the PKPD model

```{r, echo = TRUE, eval = FALSE, comment=''} 

modelEquations = list(  "Deriv_Cc" = "dose_RespPK/V*ka*exp(-ka*t) - Cl/V*Cc",
                        "Deriv_E" = "Rin*(1-Imax*(Cc**gamma)/(Cc**gamma + IC50**gamma))-kout*E" ) 
```


### Set mu and omega for each parameter

```{r, echo = TRUE, eval = FALSE, comment=''} 

modelParameters = list(
  
  ModelParameter( name = "V",
                  distribution = LogNormal( mu = 0.74, omega = 0.316 ) ),
  
  ModelParameter( name = "Cl",
                  distribution = LogNormal( mu = 0.28, omega = 0.456 ) ),
  
  ModelParameter( name = "ka",
                  distribution = LogNormal( mu = 10, omega = sqrt( 0 ) ), fixedMu = TRUE ),
  
  ModelParameter( name = "kout",
                  distribution = LogNormal( mu = 6.14, omega = 0.947 ) ),
  
  ModelParameter( name = "Rin",
                  distribution = LogNormal( mu = 614, omega = sqrt( 0 ) ), fixedMu = TRUE ),
  
  ModelParameter( name = "Imax",
                  distribution = LogNormal( mu = 0.76, omega = 0.439 ) ),
  
  ModelParameter( name = "IC50",
                  distribution = LogNormal( mu = 9.22, omega = 0.452 ) ),
  
  ModelParameter( name = "gamma",
                  distribution = LogNormal( mu = 2.77, omega = 1.761 ) ) )

```


### Create the error model to the responses PK and PD. 

```{r, echo = TRUE, eval = FALSE, comment=''} 

errorModelRespPK = Combined1( output = "RespPK", sigmaInter = 0, sigmaSlope = 0.21 )

errorModelRespPD = Constant( output = "RespPD", sigmaInter = 9.6 )

modelError = list( errorModelRespPK, errorModelRespPD )
```

### Define the arms and for each arm 
+ create the administration parameters for the response PK
+ create the sampling times for the responses PK and PD

```{r, echo = TRUE, eval = FALSE, comment=''} 


# sampling times
samplingTimesRespPK = SamplingTimes( outcome = "RespPK",
                                     samplings = c( 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4 ) )

samplingTimesRespPD = SamplingTimes( outcome = "RespPD",
                                     samplings = c( 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4 ) )

# Define the arms and the administrations

administrationRespPK1 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 0.2 ) )

arm1 = Arm( name = "0.2mg Arm",
            size = 6,
            administrations = list( administrationRespPK1 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )

administrationRespPK2 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 0.64 ) )

arm2 = Arm( name = "0.64mg Arm",
            size = 6,
            administrations = list( administrationRespPK2 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )

administrationRespPK3 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 2 ) )

arm3 = Arm( name = "2mg Arm",
            size = 6,
            administrations = list( administrationRespPK3 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )

administrationRespPK4 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 6.24 ) )

arm4 = Arm( name = "6.24mg Arm",
            size = 6,
            administrations = list( administrationRespPK4 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )

administrationRespPK5 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 11.24 ) )

arm5 = Arm( name = "11.24mg Arm",
            size = 6,
            administrations = list( administrationRespPK5 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )

administrationRespPK6 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 20 ) )

arm6 = Arm( name = "20mg Arm",
            size = 6,
            administrations = list( administrationRespPK6 ) ,
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            initialCondition = list( "Cc" = 0,
                                     "E" = 100 ) )
```

### Add the arms to the design `design1`

```{r, echo = TRUE, eval = FALSE, comment=''} 

design1 = Design( name = "design1", arms = list( arm1, arm2, arm3, arm4, arm5, arm6 ) )
```

### Evaluate the population FIM

```{r, echo = TRUE, eval = FALSE, comment=''} 

evaluationPop = Evaluation( name = " ",
                            modelEquations = modelEquations,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                            designs = list( design1 ),
                            fimType = "population",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationPop = run( evaluationPop )
```

### Display the results of the design evaluation

```{r, echo = TRUE, eval = FALSE, comment=''} 

show( evaluationPop )
fisherMatrix = getFisherMatrix( evaluationPop )
getCorrelationMatrix( evaluationPop )
getSE( evaluationPop )
getRSE( evaluationPop )
getShrinkage( evaluationPop )
getDeterminant( evaluationPop )
getDcriterion( evaluationPop )
```

### Graphs of the results of the design evaluation

```{r, echo = TRUE, eval = FALSE, comment=''} 

plotOptions = list( unitTime = c("hour"), unitOutcomes = c("mcg/mL","DI%")  )
plotEvaluation = plotEvaluation( evaluationPop, plotOptions )
plotSensitivityIndices = plotSensitivityIndices( evaluationPop, plotOptions )

SE = plotSE( evaluationPop )
RSE = plotRSE( evaluationPop )

# examples
plotOutcomesEvaluationRespPK = plotEvaluation[["design1"]][["20mg Arm"]][["RespPK"]]
plotOutcomesEvaluationRespPD = plotEvaluation[["design1"]][["20mg Arm"]][["RespPD"]]
plotSensitivityIndice_RespPK_Cl = plotSensitivityIndices[["design1"]][["20mg Arm"]][["RespPK"]][["Cl"]]
```

### Report for the design evaluation

```{r, echo = TRUE, eval = FALSE, comment=''} 

outputPath = getwd()
outputFile = "Example01_EvaluationPopFIM.html"
plotOptions = list( unitTime=c("hour"), unitOutcomes = c("mcg/mL","DI%") )
Report( evaluationPop, outputPath, outputFile, plotOptions )
```

### Evaluate the individual and Bayesian FIMs

```{r, echo = TRUE, eval = FALSE, comment=''} 

evaluationInd = Evaluation( name = " ",
                            modelEquations = modelEquations,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                            designs = list( design1 ),
                            fimType = "individual",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationInd = run( evaluationInd )

evaluationBay = Evaluation( name = " ",
                            modelEquations = modelEquations,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                            designs = list( design1 ),
                            fimType = "Bayesian",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationBay = run( evaluationBay )
```

### Display the results of the design evaluation

```{r, echo = TRUE, eval = FALSE, comment=''} 


show( evaluationInd )
show( evaluationBay )

fisherMatrix = getFisherMatrix( evaluationInd )
getCorrelationMatrix( evaluationInd )
getSE( evaluationInd )
getRSE( evaluationInd )
getShrinkage( evaluationInd )
getDeterminant( evaluationInd )
getDcriterion( evaluationInd )

fisherMatrix = getFisherMatrix( evaluationBay )

getCorrelationMatrix( evaluationBay )
getSE( evaluationBay )
getRSE( evaluationBay )
getShrinkage( evaluationBay )
getDeterminant( evaluationBay )
getDcriterion( evaluationBay )
```

### Reports of the results of the design evaluation

```{r, echo = TRUE, eval = FALSE, comment=''} 


outputPath = getwd()

plotOptions = list( unitTime=c("hour"), unitOutcomes= c("mcg/mL","DI%") )

outputFile = "Example01_EvaluationIndFIM.html"
Report( evaluationInd, outputPath, outputFile, plotOptions )

outputFile = "Example01_EvaluationBayFIM.html"
Report( evaluationBay, outputPath, outputFile, plotOptions )

```

# Design optimization

### Create an administration for response PK and samplings times

```{r, echo = TRUE, eval = FALSE, comment=''} 

administrationRespPK = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 6.24 ) )

```
We create sampling times that will be used in the initial design for comparison during the optimization process. 

```{r, echo = TRUE, eval = FALSE, comment=''} 

samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 0.25, 0.75, 1, 1.5, 2, 4, 6 ) )

samplingTimesRespPD = SamplingTimes( outcome = "RespPD", samplings = c( 0.25, 0.75, 1.5, 2, 3, 6, 8, 12 ) )
```

### Create a sampling constraint for the responses PK and PD

+ Set the vector of allowed sampling times for the responses PK and PD
+ Fix certain time values
+ Set the number of optimizable sampling times 

```{r, echo = TRUE, eval = FALSE, comment=''} 

samplingConstraintsRespPK  = SamplingTimeConstraints( outcome = "RespPK",
                                                      initialSamplings = c( 0.25, 0.75, 1, 1.5, 2, 4, 6 ),
                                                      fixedTimes = c( 0.25, 4 ),
                                                      numberOfsamplingsOptimisable = 4 )

samplingConstraintsRespPD  = SamplingTimeConstraints( outcome = "RespPD",
                                                      initialSamplings = c( 0.25, 0.75, 1.5, 2, 3, 6, 8, 12 ),
                                                      fixedTimes = c( 2, 6 ),
                                                      numberOfsamplingsOptimisable = 4 )
```
### Set the initial vectors for the FedorovWynn algorithm

```{r, echo = TRUE, eval = FALSE, comment=''} 

initialElementaryProtocols = list( c( 0.25, 0.75, 1, 4 ), c( 1.5, 2, 6, 12 ) )
```

### Create an administration constraint with the vector of allowed amount of doses for the response PK

```{r, echo = TRUE, eval = FALSE, comment=''} 

administrationConstraintsRespK = AdministrationConstraints( outcome = "RespPK", 
                                                            doses = list( 0.2, 0.64, 2, 6.24, 11.24, 20 ) )

```

### Create the constraint design to the project
+ total number of individuals to be considered
+ administration constraint
+ sampling constraint

```{r, echo = TRUE, eval = FALSE, comment=''} 


armConstraint = Arm( name = "armConstraint",
                     size = 30,
                     administrations = list( administrationRespPK ),
                     samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
                     administrationsConstraints = list( administrationConstraintsRespK ),
                     samplingTimesConstraints = list( samplingConstraintsRespPK, samplingConstraintsRespPD ),
                     initialCondition = list( "Cc" = 0,
                                              "E" = "Rin/kout"  ) )

designConstraint = Design( name = "designConstraint", 
                           arms = list( armConstraint ), numberOfArms = 30 )
```

### Set the vector of initial proportions or numbers of subjects for each elementary design

```{r, echo = TRUE, eval = FALSE, comment=''} 

numberOfSubjects = c(30)
proportionsOfSubjects = c(30)/30
```

### Set the the parameters of the FedorovWynn algorithm and run the algorithm for the design optimization with a population FIM

```{r, echo = TRUE, eval = FALSE, comment=''} 

optimizationFWPopFIM = Optimization(  name = "PKPD_ODE_multi_doses_populationFIM",
                                      modelEquations = modelEquations,
                                      modelParameters = modelParameters,
                                      modelError = modelError,
                                      
                                      optimizer = "FedorovWynnAlgorithm",
                                      
                                      optimizerParameters = list( elementaryProtocols = initialElementaryProtocols,
                                                                  numberOfSubjects = numberOfSubjects,
                                                                  proportionsOfSubjects = proportionsOfSubjects,
                                                                  showProcess = T ),
                                      
                                      designs = list( designConstraint ),
                                      
                                      fimType = "population",
                                      
                                      outputs = list( "RespPK" = "Cc","RespPD" = "E" ),
                                      
                                      odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

optimizationFWPopFIM = run( optimizationFWPopFIM )

saveRDS( optimizationFWPopFIM, "optimizationFWPopFIM.RDS" )
```

### Display and plot the results of the design optimization

```{r, echo = TRUE, eval = FALSE, comment=''} 

show( optimizationFWPopFIM )
fisherMatrix = getFisherMatrix( optimizationFWPopFIM )
getCorrelationMatrix( optimizationFWPopFIM )
getSE( optimizationFWPopFIM )
getRSE( optimizationFWPopFIM )
getShrinkage( optimizationFWPopFIM )
getDeterminant( optimizationFWPopFIM )
getDcriterion( optimizationFWPopFIM )

# plot the frequencies 
plotFrequencies = plotFrequencies( optimizationFWPopFIM )
plotFrequencies
```

### Reports for the design optimization

```{r, echo = TRUE, eval = FALSE, comment=''} 

outputFile = "Example01_OptimizationFWPopFIM.html"
Report( optimizationFWPopFIM, outputPath, outputFile, plotOptions )
```

### Set the the parameters of the Multiplicative Algorithm algorithm and run the algorithm for the design optimization with a population FIM

```{r, echo = TRUE, eval = FALSE, comment=''} 

optimizationMultPopFIM = Optimization( name = "PKPD_ODE_multi_doses_populationFIM",
                                       modelEquations = modelEquations,
                                       modelParameters = modelParameters,
                                       modelError = modelError,
                                       
                                       optimizer = "MultiplicativeAlgorithm",
                                       
                                       optimizerParameters = list( lambda = 0.99,
                                                                   numberOfIterations = 1000,
                                                                   weightThreshold = 0.01,
                                                                   delta = 1e-04, showProcess = T ),
                                       
                                       designs = list( designConstraint ),
                                       
                                       fimType = "population",
                                       
                                       outputs = list( "RespPK" = "Cc","RespPD" = "E" ),
                                       
                                       odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

```
 
### Run the Multiplicative Algorithm algorithm for the optimization with a population FIM

```{r, echo = TRUE, eval = FALSE, comment=''} 

optimizationMultPopFIM = run( optimizationMultPopFIM )
saveRDS( optimizationMultPopFIM, "optimizationMultPopFIM.RDS" )

```

### Display and plot the results of the design optimization

```{r, echo = TRUE, eval = FALSE, comment=''} 

show( optimizationMultPopFIM )
fisherMatrix = getFisherMatrix( optimizationMultPopFIM )
getCorrelationMatrix( optimizationMultPopFIM )
getSE( optimizationMultPopFIM )
getRSE( optimizationMultPopFIM )
getShrinkage( optimizationMultPopFIM )
getDeterminant( optimizationMultPopFIM )
getDcriterion( optimizationMultPopFIM )

# plot the weight
plotWeights = plotWeights( optimizationMultPopFIM )
plotWeights

```

### Create and save the report for the design optimization

```{r, echo = TRUE, eval = FALSE, comment=''} 

outputPath = getwd()
plotOptions = list( unitTime = c( "hour" ), unitOutcomes= c( "mcg/mL" , "DI%" ) )
outputFile = "Example01_OptimizationMultPopFIM.html"
Report( optimizationMultPopFIM, outputPath, outputFile, plotOptions )
```

# References  
```{r global_options_end, echo = FALSE, include = FALSE}
options(backup_options)
```