Use Case 1 Optimal design evaluation

This example script is intended to to illustrate optimal design with PFIM and Pop ED

The following steps are implemented in this workflow:

Prepared by Mike K Smith, Pfizer Ltd, Sandwich, UK 04 August 2016

Initialisation

Clear workspace and set working directory

rm(list=ls(all=F))
getwd()
## [1] "C:/SEE/MDL_IDE/workspace/UseCasesDemo/scripts/Design"

Some useful functions for working with the SEE and workspaces in the MDL-IDE A system variable has been set which retrieves the directory of the workspaces

Sys.getenv("MDLIDE_WORKSPACE_HOME")
## [1] "C:\\SEE\\MDL_IDE\\workspace"

By typing “~” you can get to this directory

mydir <- file.path("~","UseCasesDemo","models","Design")
setwd(mydir)

Another function from the ddmore R package retrieves the directory that the SEE is installed in. You can then navigate relative to these directories.

ddmore:::DDMORE.checkConfiguration()
## [1] "C:/SEE"

Set name of .mdl file

mdlfile <- "UseCase1_design1_optFW.mdl"
model <- tools::file_path_sans_ext(mdlfile)

wd <- file.path(mydir,model)
dir.create(wd)

Copy the dataset and the .mdl file available under “models” into the working directory

file.copy(file.path(mydir,mdlfile),wd)
## [1] TRUE

Set the working directory.

setwd(wd)

The working directory needs to be set differently when knitr R package is used to spin the file

library(knitr)
opts_knit$set(root.dir = wd) 
cache.path <- file.path(wd, 'cache')
dir.create(cache.path)
opts_chunk$set(cache.path=cache.path)
figure.path <- file.path(wd, 'figure')
dir.create(figure.path)
opts_chunk$set(figure.path=figure.path)

List files available in working directory

list.files()
## [1] "cache"                      "figure"                    
## [3] "UseCase1_design1_optFW.mdl"

Read the different MDL objects needed for upcoming tasks. An MDL file may contain more than one object of any type (though typically only one model!) so these functions return a list of all objects of that type found in the target MDL file. To pick out the first, use the double square bracket notation in R to retrieve the first item of the list.

myModelObj <- getModelObjects(mdlfile)[[1]] 
myDesignObj <- getDesignObjects(mdlfile)[[1]]
myParameterObj <- getParameterObjects(mdlfile)[[1]]
myTaskPropertiesObj <- getTaskPropertiesObjects(mdlfile)[[1]]

Design evaluation using PFIM

The Design Object settings and Task Properties settings dictate whether evaluation or optimisation of the trial design is performed. The runPFIM function takes an MDL file, converts to PharmML, runs the converter and then runs a created batch script to run PFIM.

pharmmlfile <- as.PharmML(mdlfile)

The runPFIM function takes the MDL file (or converted PharmML) as inputs, and runs the converter to create the required PFIM .R files and an associated .bat file for running PFIM in R. Note that PFIM is run in an independent R instance. An option “run=FALSE” will perform the conversion and stop.

runPFIM(pharmmlfile=pharmmlfile, jarLocation=file.path(ddmore:::DDMORE.checkConfiguration(),"PFIM","converter"))
## [1] "java -jar C:/SEE/PFIM/converter/pfim.jar -p C:/SEE/PFIM/PFIM4.0/Program -i C:\\SEE\\MDL_IDE\\workspace\\UseCasesDemo\\models\\Design\\UseCase1_design1_optFW\\UseCase1_design1_optFW.xml -o C:/SEE/MDL_IDE/workspace/UseCasesDemo/models/Design/UseCase1_design1_optFW/PFIM"
myOutputFile <- myTaskPropertiesObj@OPTIMISE$blocks[[1]]$TARGET_SETTINGS$output
myOutputFile <- gsub("\\\"","",myOutputFile)
cat(readLines(file.path(getwd(),"PFIM",myOutputFile) ),sep="\n")
## PFIM 4.0  
##  
## Option: 1 
## 
##  
## Project: warfarin_PK_ODE
##  
## Date: Fri Aug 19 16:11:57 2016
##  
## 
##  
## **************************** INPUT SUMMARY ********************************
##  
## Differential Equations form of the model:  
##  
## function (t, y, p) 
## {
##     CL <- p[1]
##     V <- p[2]
##     KA <- p[3]
##     RATEIN <- ((y[1]) * (KA))
##     CC <- ((y[2])/(V))
##     yd1 <- -(RATEIN)
##     yd2 <- RATEIN - ((((CL) * (y[2])))/(V))
##     return(list(c(yd1, yd2), CC))
## }
## 
##  
## Initial Conditions at time 0 : 
##  
## 100 0 
## 
## Error tolerance for solving differential equations system: RtolEQ = 1e-08 , AtolEQ = 1e-08 , Hmax =  Inf
##  
## Initial design: 
## 
##  
## Sample times for response: A 
##                 Protocol subjects  condinit
## 1 c=(1e-04, 36, 96, 120)       33 c(100, 0)
## 
##  
## Total number of samples: 132
##  
## Associated criterion value: 57.1694
##  
## Identical sampling times for each response: TRUE
##  
## Random effect model: Trand =  2
##  
## Variance error model response A : ( 0.1 + 0.1 *f)^2
## 
##  
## 
## Optimization step:  
##  
## Sampling windows for the response: A 
## Window 1 : t= 1e-04 2 6 12 24 36 48 72 96 120 
##     Nb of sampling points to be taken in this window, n[ 1 ]= 4 5 
## Maximum total number of points in one elementary protocol : 5 
## Minimum total number of points in one elementary protocol : 4 
## 
##  
## 
## Now evaluating the Fisher Information Matrix for the 462 protocols generated 
## 
##  
## BEST ONE GROUP PROTOCOL: 
##  
## Sample times for response: A 
##                      times freq Subjects  condinit
## 1 c(1e-04, 2, 12, 24, 120)    1     26.4 c(100, 0)
## 
##  
## Associated criterion: 1117.465
##  
## 
##  
## **************************** OPTIMISED DESIGN *****************************
##  
## 
##  
## Optimised design: 
## Sample times for response: A 
##                      times      freq  Subjects  condinit
## 1        c(2, 12, 24, 120) 0.4852886 15.574730 c(100, 0)
## 2     c(1e-04, 2, 24, 120) 0.4017606 12.894002 c(100, 0)
## 3 c(1e-04, 2, 12, 24, 120) 0.1129508  3.625014 c(100, 0)
## 
##  
## 
##  
## Associated optimised criterion: 1153.623
##  
## 
##  
## FIM saved in FIM.txt
##  
##  
## Computation of the Population Fisher information matrix: option =  1
##  
## ******************* FISHER INFORMATION MATRIX ******************
##  
##              [,1]       [,2]       [,3]         [,4]         [,5]         [,6]
## [1,] 29676.979142  -3.604482  -47.55527 0.000000e+00    0.0000000   0.00000000
## [2,]    -3.604482   4.565460  -11.03317 0.000000e+00    0.0000000   0.00000000
## [3,]   -47.555267 -11.033171 1629.46077 0.000000e+00    0.0000000   0.00000000
## [4,]     0.000000   0.000000    0.00000 1.372112e+03    0.1767931   0.09485664
## [5,]     0.000000   0.000000    0.00000 1.767931e-01 1331.1086820  16.96292642
## [6,]     0.000000   0.000000    0.00000 9.485664e-02   16.9629264 711.43842953
## [7,]     0.000000   0.000000    0.00000 5.361802e+01   22.9657175  82.17412927
## [8,]     0.000000   0.000000    0.00000 1.693120e+02  202.4316982 596.14402209
##            [,7]      [,8]
## [1,]    0.00000    0.0000
## [2,]    0.00000    0.0000
## [3,]    0.00000    0.0000
## [4,]   53.61802  169.3120
## [5,]   22.96572  202.4317
## [6,]   82.17413  596.1440
## [7,] 3346.93671  416.8145
## [8,]  416.81448 3916.3378
## 
##  
## 
##  
## ************************** EXPECTED STANDARD ERRORS ************************
##  
## ------------------------ Fixed Effects Parameters -------------------------
##  
##     Beta   StdError      RSE  
## CL 0.100 0.00580531 5.805310 %
## V  8.000 0.47191688 5.898961 %
## KA 0.362 0.02497901 6.900278 %
## 
##  
## 
##  
## ------------------------- Variance of Inter-Subject Random Effects ----------------------
##  
##    omega2   StdError      RSE  
## CL    0.1 0.02708381 27.08381 %
## V     0.1 0.02752124 27.52124 %
## KA    0.1 0.04016266 40.16266 %
## 
##  
## ------------------------ Standard deviation of residual error ------------------------ 
##  
##            Sigma   StdError      RSE  
## sig.interA   0.1 0.01740503 17.40503 %
## sig.slopeA   0.1 0.01731622 17.31622 %
## 
##  
## ******************************* DETERMINANT ********************************
##  
## 3.136988e+24
##  
## ******************************** CRITERION *********************************
##  
## 1153.624
##  
## 
##  
## 
##  
## ******************* EIGENVALUES OF THE FISHER INFORMATION MATRIX ******************
##  
##         FixedEffects VarianceComponents
## min      3144.974032           4.490027
## max     29677.060207        3144.974032
## max/min     9.436345         700.435403
## 
##  
## ******************* CORRELATION MATRIX ******************
##  
##            [,1]       [,2]       [,3]         [,4]          [,5]        [,6]
## [1,] 1.00000000 0.01075584 0.00815865  0.000000000  0.0000000000  0.00000000
## [2,] 0.01075584 1.00000000 0.12799541  0.000000000  0.0000000000  0.00000000
## [3,] 0.00815865 0.12799541 1.00000000  0.000000000  0.0000000000  0.00000000
## [4,] 0.00000000 0.00000000 0.00000000  1.000000000  0.0068333634  0.02822632
## [5,] 0.00000000 0.00000000 0.00000000  0.006833363  1.0000000000  0.01549536
## [6,] 0.00000000 0.00000000 0.00000000  0.028226320  0.0154953556  1.00000000
## [7,] 0.00000000 0.00000000 0.00000000 -0.017146959 -0.0009977094 -0.01356847
## [8,] 0.00000000 0.00000000 0.00000000 -0.076310737 -0.0879770951 -0.35478348
##               [,7]        [,8]
## [1,]  0.0000000000  0.00000000
## [2,]  0.0000000000  0.00000000
## [3,]  0.0000000000  0.00000000
## [4,] -0.0171469585 -0.07631074
## [5,] -0.0009977094 -0.08797710
## [6,] -0.0135684730 -0.35478348
## [7,]  1.0000000000 -0.10094369
## [8,] -0.1009436857  1.00000000
## 
## 
##  
## Time difference of 1.975769 secs
## sys.self 
##     0.85 
## 
## 

Design evaluation using PopED

PopED function as.poped takes a PharmML file and creates a poped.db database object ready for use with PopED. We can then use PopED functions directly (natively) in R.

library(PopED)

as.poped(pharmmlfile)
## Warning: PopED Warning: No PopED operation algorithm found in design step
## at line 459
## Warning: PopED Warning: No PopED-specific settings could be retrieved from
## modelling steps at line 392

create plot of model without variability

plot_model_prediction(poped.db)
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 85.9832
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work
## (> maxsteps ) was done, but integration was not successful - increase
## maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).

plot of chunk Plot_model_predictions_from_PopED

get predictions from model

popED.pred <- model_prediction(poped.db)
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 85.9832
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work
## (> maxsteps ) was done, but integration was not successful - increase
## maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
popED.pred[order(popED.pred$Time),]
##      Time         PRED Group Model DOSE_1_AMT DOSE_1_TIME
## 1 1.0e-04 0.0004524915     1     1        100           0
## 2 3.6e+01 8.2553862095     1     1        100           0
## 3 9.6e+01           NA     1     1        100           0
## 4 1.2e+02           NA     1     1        100           0

evaluate initial design

#FIM <- POPED_OPTIM(POPED.DB, OPT_XT=T) 
#FIM
#DET(FIM)
#GET_RSE(FIM,POPED.DB)