Guided-demo: Exploring 'ddmore' R package functionality of time to event models

UseCase14_1 : ODE + time to event model for exact and right censored information

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

Initialisation

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"

Model Development

ESTIMATE model parameters using Monolix

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

Estimation with NONMEM

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.

Simulation using simulx

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)

plot of chunk unnamed-chunk-16