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_eval.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_eval.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_eval\\UseCase1_design1_eval.xml -o C:/SEE/MDL_IDE/workspace/UseCasesDemo/models/Design/UseCase1_design1_eval/PFIM"
myOutputFile <- myTaskPropertiesObj@EVALUATE$blocks[[1]]$TARGET_SETTINGS$output
myOutputFile <- gsub("\\\"","",myOutputFile)
cat(readLines(file.path(getwd(),"PFIM",myOutputFile) ),sep="\n")
## PFIM 4.0 
##  
## Project:  warfarin_PK_ODE
##  
## Date:  Fri Aug 19 16:03:22 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))
## }
## 
##  
## Design:  
## Sample times for response: A 
##                               times subjects
## 1 c(1e-04, 24, 36, 48, 72, 96, 120)       33
## 
##  
## Initial Conditions at time 0 : 
##  
## 100 0 
## 
##  
## Random effect model: Trand =  2
##  
## Variance error model response A : ( 0.1 + 0.1 *f)^2
## 
##  
## Error tolerance for solving differential equations system: RtolEQ = 1e-08 , AtolEQ = 1e-08 , Hmax =  Inf
##  
## Computation of the Population Fisher information matrix: option =  1
##  
## ******************* FISHER INFORMATION MATRIX ******************
##  
##              [,1]     [,2]       [,3]         [,4]         [,5]        [,6]
## [1,] 31893.526555 1.086082 315.075120 0.000000e+00 0.000000e+00 0.000000000
## [2,]     1.086082 4.685126   3.590723 0.000000e+00 0.000000e+00 0.000000000
## [3,]   315.075120 3.590723   5.897095 0.000000e+00 0.000000e+00 0.000000000
## [4,]     0.000000 0.000000   0.000000 1.541208e+03 1.143829e-02 1.971067179
## [5,]     0.000000 0.000000   0.000000 1.143829e-02 1.362256e+03 1.638388900
## [6,]     0.000000 0.000000   0.000000 1.971067e+00 1.638389e+00 0.009048292
## [7,]     0.000000 0.000000   0.000000 2.159723e+01 3.514215e+01 0.064850000
## [8,]     0.000000 0.000000   0.000000 8.137515e+01 2.355279e+02 0.425930721
##            [,7]         [,8]
## [1,]    0.00000 0.000000e+00
## [2,]    0.00000 0.000000e+00
## [3,]    0.00000 0.000000e+00
## [4,]   21.59723 8.137515e+01
## [5,]   35.14215 2.355279e+02
## [6,]    0.06485 4.259307e-01
## [7,] 7272.98302 3.375624e+03
## [8,] 3375.62444 1.903638e+04
## 
##  
## ************************** EXPECTED STANDARD ERRORS ************************
##  
## ------------------------ Fixed Effects Parameters -------------------------
##  
##     Beta   StdError        RSE  
## CL 0.100 0.04489998   44.89998 %
## V  8.000 3.48569370   43.57117 %
## KA 0.362 4.52143436 1249.01502 %
## 
##  
## ------------------------- Variance of Inter-Subject Random Effects ----------------------
##  
##    omega2    StdError         RSE  
## CL    0.1  0.03174683    31.74683 %
## V     0.1  0.03244822    32.44822 %
## KA    0.1 14.81373634 14813.73634 %
## 
##  
## ------------------------ Standard deviation of residual error ------------------------ 
##  
##            Sigma    StdError       RSE  
## sig.interA   0.1 0.012240439 12.240439 %
## sig.slopeA   0.1 0.007574304  7.574304 %
## 
##  
## ******************************* DETERMINANT ********************************
##  
## 8.863853e+15
##  
## ******************************** CRITERION *********************************
##  
## 98.50376
##  
## 
##  
## ******************* EIGENVALUES OF THE FISHER INFORMATION MATRIX ******************
##  
##         FixedEffects VarianceComponents
## min      6373.288687       4.556899e-03
## max     31896.639484       1.540831e+03
## max/min     5.004738       3.381315e+05
## 
##  
## ******************* CORRELATION MATRIX ******************
##  
##            [,1]       [,2]       [,3]         [,4]          [,5]         [,6]
## [1,]  1.0000000  0.9833930 -0.9921931  0.000000000  0.0000000000  0.000000000
## [2,]  0.9833930  1.0000000 -0.9911775  0.000000000  0.0000000000  0.000000000
## [3,] -0.9921931 -0.9911775  1.0000000  0.000000000  0.0000000000  0.000000000
## [4,]  0.0000000  0.0000000  0.0000000  1.000000000  0.3279993284 -0.596716225
## [5,]  0.0000000  0.0000000  0.0000000  0.327999328  1.0000000000 -0.548906304
## [6,]  0.0000000  0.0000000  0.0000000 -0.596716225 -0.5489063037  1.000000000
## [7,]  0.0000000  0.0000000  0.0000000 -0.003077763  0.0006540826  0.002154311
## [8,]  0.0000000  0.0000000  0.0000000 -0.008307933 -0.0350477964 -0.004591849
##               [,7]         [,8]
## [1,]  0.0000000000  0.000000000
## [2,]  0.0000000000  0.000000000
## [3,]  0.0000000000  0.000000000
## [4,] -0.0030777633 -0.008307933
## [5,]  0.0006540826 -0.035047796
## [6,]  0.0021543106 -0.004591849
## [7,]  1.0000000000 -0.286639404
## [8,] -0.2866394039  1.000000000
## 
## 
##  
## Time difference of 0.7308791 secs

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 426
## Warning: PopED Warning: No PopED-specific settings could be retrieved from
## modelling steps at line 359

create plot of model without variability

plot_model_prediction(poped.db)

plot of chunk Plot_model_predictions_from_PopED

get predictions from model

popED.pred <- model_prediction(poped.db)

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 2.4e+01 9.5892404798     1     1        100           0
## 3 3.6e+01 8.2553862095     1     1        100           0
## 4 4.8e+01 7.1055007860     1     1        100           0
## 5 7.2e+01 5.2638847219     1     1        100           0
## 6 9.6e+01 3.8995817136     1     1        100           0
## 7 1.2e+02 2.8888811865     1     1        100           0

evaluate initial design

FIM <- evaluate.fim(poped.db) 
FIM
##              [,1]        [,2]        [,3]         [,4]         [,5]
## [1,] 7965.4235899  3.12636864 80.90450214 -467.7882762 -99.35817267
## [2,]    3.1263686  0.82528370  0.65624709   -2.8324658  -0.30263736
## [3,]   80.9045021  0.65624709  1.31214135   -6.9495290  -1.17690220
## [4,] -467.7882762 -2.83246578 -6.94952897  621.1957316 109.39784174
## [5,]  -99.3581727 -0.30263736 -1.17690220  109.3978417 222.15111005
## [6,]   -0.6727683 -0.01069187 -0.01124996    0.4624558   0.06474483
## [7,] -247.8071581 -2.77145103 -7.14362482  217.3891712  32.04991219
## [8,] -257.5734156 -5.66111959 -1.48402595  148.1571304  38.74891733
##             [,6]        [,7]        [,8]
## [1,] -0.67276829 -247.807158 -257.573416
## [2,] -0.01069187   -2.771451   -5.661120
## [3,] -0.01124996   -7.143625   -1.484026
## [4,]  0.46245580  217.389171  148.157130
## [5,]  0.06474483   32.049912   38.748917
## [6,]  0.01268646    1.600999    3.714749
## [7,]  1.60099870 1257.992222    1.190294
## [8,]  3.71474943    1.190294 2806.688130
det(FIM)
## [1] 21921978755
get_rse(FIM,poped.db)
##      bpop[1]      bpop[2]      bpop[3]      bpop[4]      bpop[5] 
##    252.55288    244.91357   7072.59358     44.42654     70.33416 
##       D[1,1]       D[2,2]       D[3,3] 
## 132663.42851    610.58447    569.98209