This example script is intended to illustrate how to use the 'ddmore' R package to perform a M&S workflow using the DDMoRe Standalone Execution Environment (SEE).
The following steps are implemented in this workflow:
To run a task, select with the cursor any code lines you wish to execute and press CTRL+R+R in your keyboard. An HTML file containing the commands in this file and associated output will be provided to allow the user to compare the results
Clear workspace and set working directory under 'UsesCasesDemo' project
rm(list=ls(all=FALSE))
mydir <- file.path(Sys.getenv("MDLIDE_WORKSPACE_HOME"),"UseCasesDemo")
setwd(mydir)
Create a working directory under 'models' folder where results are stored
uc<-"UseCase14_1"
datafile <- "warfarin_TTE_exact.csv"
mdlfile <- paste0(uc,".mdl")
Create a new folder to be used as working directory and where results will be stored
wd <- file.path(mydir,uc)
dir.create(wd)
Copy the dataset and the .mdl file available under “models” into the working directory
file.copy(file.path(mydir,"models", datafile),wd)
## [1] TRUE
file.copy(file.path(mydir,"models",mdlfile),wd)
## [1] TRUE
Set the working directory.
setwd(file.path(mydir,uc))
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 = file.path(mydir,uc))
List files available in working directory
list.files()
## [1] "UseCase14_1.mdl" "warfarin_TTE_exact.csv"
The ddmore “estimate” function translates the contents of the .mdl file to a target language and then estimates parameters using the target software. After estimation, the output is converted to a Standardised Output object which is saved in a .SO.xml file.
Translated files and Monolix output will be returned in the ./Monolix subfolder. The Standardised Output object (.SO.xml) is read and parsed into an R object called “mlx” of (S4) class “StandardOutputObject”.
mlx <- estimate(mdlfile, target="MONOLIX", subfolder="Monolix")
## -- Tue Aug 16 18:10:35 2016
## New
## Submitted
## Job bd8c81f2-e002-4cd6-b793-325683397911 progress:
## Running [ ............................... ]
## Importing Results
## Copying the result data back to the local machine for job ID bd8c81f2-e002-4cd6-b793-325683397911...
## From C:\Users\zparra\AppData\Local\Temp\Rtmpk5tG0t\DDMORE.job1c847de06383 to D:/SEE-Prod5_RC4/MDL_IDE/workspace/UseCasesDemo/UseCase14_1/Monolix
## Done.
##
##
## The following main elements were parsed successfully:
## ToolSettings
## RawResults
## Estimation::PopulationEstimates::MLE
## Estimation::PrecisionPopulationEstimates::MLE
## Estimation::IndividualEstimates::Estimates
## Estimation::IndividualEstimates::RandomEffects
## Estimation::OFMeasures::IndividualContribToLL
## Estimation::OFMeasures::InformationCriteria
## Estimation::OFMeasures::LogLikelihood
##
## Completed
## -- Tue Aug 16 18:21:01 2016
slotNames(mlx)
## [1] "ToolSettings" "RawResults" "TaskInformation"
## [4] "Estimation" "ModelDiagnostic" "Simulation"
## [7] "OptimalDesign" ".pathToSourceXML"
The ddmore “LoadSOObj” reads and parsed existing Standardise Output Objects
mlx <- LoadSOObject(“Monolix/UseCase14_1.SO.xml”)
The ddmore function “getPopulationParameters” extracts the Population Parameter values from an R object of (S4) class “StandardOutputObject” and returns the estimates. See documentation for getPopulationParameters to see other arguments and settings for this function.
parameters_mlx <- getPopulationParameters(mlx, what="estimates")$MLE
print(parameters_mlx)
## POP_BTATRT POP_HBASE
## 0.03009 9.88372
print(getPopulationParameters(mlx, what="precisions"))
## $MLE
## Parameter MLE SE RSE
## 1 POP_BTATRT 0.03009 0.10203 339.07
## 2 POP_HBASE 9.88372 2.14151 21.67
print(getEstimationInfo(mlx))
## $OFMeasures
## $OFMeasures$LogLikelihood
## $OFMeasures$LogLikelihood[[1]]
## [1] -255.165
##
##
## $OFMeasures$IndividualContribToLL
## Subject ICtoLL
## 1 1 -0.76
## 2 2 -0.78
## 3 3 -3.69
## 4 4 -3.89
## 5 5 -3.64
## 6 6 -3.58
## 7 7 -0.80
## 8 8 -3.88
## 9 9 -0.76
## 10 10 -4.18
## 11 11 -4.22
## 12 12 -3.98
## 13 13 -0.85
## 14 14 -3.53
## 15 15 -4.30
## 16 16 -0.76
## 17 17 -0.80
## 18 18 -4.04
## 19 19 -4.17
## 20 20 -4.19
## 21 21 -0.80
## 22 22 -4.17
## 23 23 -3.74
## 24 24 -0.78
## 25 25 -0.78
## 26 26 -0.85
## 27 27 -0.80
## 28 28 -3.80
## 29 29 -3.59
## 30 30 -0.76
## 31 31 -4.28
## 32 32 -0.85
## 33 33 -0.78
## 34 34 -3.68
## 35 35 -0.80
## 36 36 -4.05
## 37 37 -4.28
## 38 38 -4.19
## 39 39 -3.69
## 40 40 -3.98
## 41 41 -3.66
## 42 42 -0.85
## 43 43 -0.80
## 44 44 -3.66
## 45 45 -3.69
## 46 46 -3.92
## 47 47 -0.85
## 48 48 -0.85
## 49 49 -0.78
## 50 50 -4.33
## 51 51 -3.67
## 52 52 -0.76
## 53 53 -0.85
## 54 54 -0.78
## 55 55 -0.76
## 56 56 -4.11
## 57 57 -3.89
## 58 58 -3.93
## 59 59 -0.76
## 60 60 -4.12
## 61 61 -3.61
## 62 62 -0.78
## 63 63 -0.78
## 64 64 -3.61
## 65 65 -0.85
## 66 66 -3.77
## 67 67 -0.78
## 68 68 -0.78
## 69 69 -3.77
## 70 70 -3.95
## 71 71 -0.78
## 72 72 -3.61
## 73 73 -3.95
## 74 74 -3.64
## 75 75 -3.91
## 76 76 -0.78
## 77 77 -4.11
## 78 78 -0.76
## 79 79 -4.15
## 80 80 -0.78
## 81 81 -4.19
## 82 82 -4.02
## 83 83 -0.78
## 84 84 -0.80
## 85 85 -0.78
## 86 86 -0.76
## 87 87 -0.76
## 88 88 -4.21
## 89 89 -4.26
## 90 90 -4.07
## 91 91 -3.74
## 92 92 -3.61
## 93 93 -0.85
## 94 94 -4.04
## 95 95 -0.85
## 96 96 -3.74
## 97 97 -4.18
## 98 98 -4.27
## 99 99 -0.76
## 100 100 -0.85
##
## $OFMeasures$InformationCriteria
## $OFMeasures$InformationCriteria$AIC
## [1] 514.33
##
## $OFMeasures$InformationCriteria$BIC
## [1] 519.54
##
##
##
## $Messages
## list()
NONMEM estimation is not possible for this use case. In the current version of the framework, there is no place to define initial conditions for an ODE which are not structural parameters so that the model runs interoperably. See UC14_2 for an alternative implementation that can be used with NONMEM.
The mlxR package has been developed to visualize and explore models that are encoded in MLXTRAN or PharmML. The ddmore function as.PharmML translates an MDL file (extension .mdl) to its PharmML representation. The output file (extension .xml) is saved in the working directory.
myPharmML <- as.PharmML(mdlfile)
Use parameter values from the Monolix estimation
parValues <- getPopulationParameters(mlx, what="estimates")$MLE
parValues <- parValues[c("POP_BTATRT","POP_HBASE")]
Add the treatment covariate
p1 <- c(parValues, TRT = 1)
p2 <- c(parValues, TRT = 4)
Parameter values used in simulation
print(p1)
## POP_BTATRT POP_HBASE TRT
## 0.03009 9.88372 1.00000
print(p2)
## POP_BTATRT POP_HBASE TRT
## 0.03009 9.88372 4.00000
Simulate events
h <- list(name='HAZ', time=seq(0, 60, by=1))
e <- list(name='Y', time=0)
Simulate 100 subjects within each treatment group
g1 <- list( size = 100, level = 'longitudinal', parameter = p1)
g2 <- list( size = 100, level = 'longitudinal', parameter = p2)
Call simulx
res <- simulx(model = myPharmML,
group = list(g1,g2),
output = list(h,e))
Kaplan-Meyer plot of the simulated time-to-event data
pl1 <- kmplotmlx(res$Y, level=0.9)
print(pl1)