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
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]]
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
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)
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