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_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]]
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
##
##
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).
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)