Installation Test

This example script is intended to to test installation of software and tools using the DDMoRe Standalone Execution Environment (SEE).

The following steps are implemented in this workflow:

Initialisation

Clear workspace and set working directory

rm(list=ls(all=F))
getwd()
## [1] "C:/SEE/MDL_IDE/workspace/InstallationTest/models"

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("~/InstallationTest/models")
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 <- "InstallationTest_DEQ.mdl"
model <- tools::file_path_sans_ext(mdlfile)

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]] 
myDataObj <- getDataObjects(mdlfile)[[1]]
myParameterObj <- getParameterObjects(mdlfile)[[1]]
myPriorObj <- getPriorObjects(mdlfile)[[1]]
myDesignObj <- getDesignObjects(mdlfile)[[1]]

You can also refer to objects by their name.

myTaskPropertiesObj_saem <- getTaskPropertiesObjects(mdlfile)[["SAEM_task"]]
myTaskPropertiesObj_focei <- getTaskPropertiesObjects(mdlfile)[["NONMEM_task"]]
myTaskPropertiesObj_mcmc <- getTaskPropertiesObjects(mdlfile)[["BUGS_task"]]
myTaskPropertiesObj_Evaluate<- getTaskPropertiesObjects(mdlfile)[["EVALUATE_task"]]
#myTaskPropertiesObj_Optimise <- getTaskPropertiesObjects(mdlfile)[["warfarin_Optimise_Explicit"]]

Exploratory Data Analysis

getDataObjects reads and parses the MDL. It doesn't actually read the data file. To do so, use the ddmore function readDataObj() to create an R object (data frame) based on the MDL data object. This uses data column headers defined in the Data Object.

myData <- readDataObj(myDataObj)

Let's look at the first 6 lines of the data set

head(myData)
##   ID AMT TIME    DV MDV
## 1  1 100    0    NA   1
## 2  1  NA    1 9.436   0
## 3  1  NA    2 8.214   0
## 4  1  NA    4 7.945   0
## 5  1  NA    8 4.978   0
## 6  1  NA   12 3.853   0
myEDAData <- myData[myData$MDV==0,]

library(ggplot2)

Plot the data using ggplot2 library

plot1 <- ggplot(aes(y=DV,x=TIME, group=ID), data=myEDAData) +
                geom_line()+
                geom_point() + 
                labs(y="Concentration" ,x="Time (hours)")
print(plot1)

plot of chunk ExploratoryAnalysis

Model Development

ESTIMATE model parameters using Monolix

You can pass MDL files (or PharmML files) into ddmore functions to perform tasks. But you can also define new Model Object Groups (MOGs) and then write these to file. This allows you to combine different elements from the MOG to do specific tasks.

To define the new MOG for estimation in Monolix:

MLX.MOG <- createMogObj(dataObj = myDataObj, 
        parObj = myParameterObj, 
        mdlObj = myModelObj, 
        taskObj = myTaskPropertiesObj_saem)

We can then write the MOG back out to an .mdl file.

mdlfile.MLX <- paste0(model,"_MLX.mdl")
writeMogObj(MLX.MOG,mdlfile.MLX)

The ddmore “estimate” function translates the contents of the .mdl file to MLXTRAN and then estimates parameters using Monolix. After estimation, the output from Monolix 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.MLX, target="MONOLIX", subfolder="Monolix")
## -- Mon Aug 15 13:56:53 2016
## New
## Submitted
## Job 2e26af2c-f6a1-4f4d-9c3b-1dbca160c823 progress:
## Running [ ... ]
## Importing Results
## Copying the result data back to the local machine for job ID 2e26af2c-f6a1-4f4d-9c3b-1dbca160c823...
## From C:\Users\smith_mk\AppData\Local\Temp\4\RtmpYF4SFO\DDMORE.job35644c3f557e to C:/SEE/MDL_IDE/workspace/InstallationTest/models/Monolix
## Done.
## 
## 
## The following main elements were parsed successfully:
##   ToolSettings
##   RawResults
##   Estimation::PopulationEstimates::MLE
##   Estimation::PrecisionPopulationEstimates::MLE
##   Estimation::IndividualEstimates::Estimates
##   Estimation::IndividualEstimates::RandomEffects
##   Estimation::Residuals::ResidualTable
##   Estimation::Predictions
##   Estimation::OFMeasures::IndividualContribToLL
##   Estimation::OFMeasures::InformationCriteria
##   Estimation::OFMeasures::LogLikelihood
## 
## Completed
## -- Mon Aug 15 13:57:56 2016
slotNames(MLX)
## [1] "ToolSettings"     "RawResults"       "TaskInformation" 
## [4] "Estimation"       "ModelDiagnostic"  "Simulation"      
## [7] "OptimalDesign"    ".pathToSourceXML"

Retrieve estimated parameters

The ddmore function getPopulationParameters extracts the Population Parameter values from an R object of (S4) class “StandardOutputObject” and returns the MLE estimates. See documentation for getPopulationParameters to see other arguments and settings for this function.

getPopulationParameters(MLX)$MLE
## Warning in getMLEPopulationParameters(SOObject, what = what): Tried to fetch the parameter interval values, however section Estimation::PrecisionPopulationEstimates::MLE::AsymptoticCI was not found in the SO Object.
##  Omitting interval values for MLE section in returned output.
##   Parameter     MLE      SE   RSE
## 1  OMEGA_CL 0.01359 0.00686 50.51
## 2   OMEGA_V 0.01232 0.00545 44.22
## 3    POP_CL 0.81344 0.03042  3.74
## 4     POP_V 8.07709 0.26952  3.34
## 5    SD_ADD 0.41413 0.04224 10.20
parameters_mlx <- getPopulationParameters(MLX, what="estimates")$MLE
parameters_mlx
##    POP_V   POP_CL  OMEGA_V OMEGA_CL   SD_ADD 
##  8.07709  0.81344  0.01232  0.01359  0.41413

Other functions exist to extract information from the SO, but the SO is also available as an R object and so it is relatively straightforward to retrieve elements directly using R commands.

Perform basic model diagnostics

Use ddmore function as.xpdb() to create an Xpose database object from an R object of (S4) class “StandardOutputObject”. We can then call Xpose functions referencing this mlx.xpdb object as the input.

datafile <- myDataObj@SOURCE$srcfile$file
mlx.xpdb<-as.xpdb(MLX,datafile)
## 
## Removed dose rows in rawData slot of SO to enable merge with Predictions data.
print(basic.gof(mlx.xpdb))

plot of chunk Basic_Goodness_of_Fit_following_Monolix_Estimation

mlx.xpdb@Data$ID <- as.numeric(as.character(mlx.xpdb@Data$ID))
print(ind.plots(mlx.xpdb))

plot of chunk Individual_Plots_following_Monolix_Estimation

ESTIMATE model parameters in NONMEM

Assemble the new Modelling Object Group (MOG) using Nonmem task properties. Note that the only change is in the Task Properties Object. All other MDL objects remain the same.

NM.MOG <- createMogObj(dataObj = myDataObj, 
        parObj = myParameterObj, 
        mdlObj = myModelObj, 
        taskObj = myTaskPropertiesObj_focei)

We can then write the MOG back out to an .mdl file as before.

mdlfile.NM <- paste0(model,"_NM.mdl")
writeMogObj(NM.MOG,mdlfile.NM)

The ddmore “estimate” function translates the contents of the .mdl file to NMTRAN and then estimates parameters using NONMEM. After estimation, the output from NONMEM is converted to a Standardised Output object which is saved in a .SO.xml file.

Translated files and NONMEM output will be returned in the ./NONMEM subfolder. The Standardised Output object (.SO.xml) is read and parsed into an R object called “nm” of (S4) class “StandardOutputObject”.

NM <- estimate(mdlfile.NM, target="NONMEM", subfolder="NONMEM")
## -- Mon Aug 15 13:57:59 2016
## New
## Submitted
## Job b7b3fbfb-d1d9-492d-8b8c-f3b8fd83f11a progress:
## Running [ .. ]
## Importing Results
## Copying the result data back to the local machine for job ID b7b3fbfb-d1d9-492d-8b8c-f3b8fd83f11a...
## From C:\Users\smith_mk\AppData\Local\Temp\4\RtmpYF4SFO\DDMORE.job3564561c1615 to C:/SEE/MDL_IDE/workspace/InstallationTest/models/NONMEM
## Done.
## 
## 
## The following main elements were parsed successfully:
##   RawResults
##   Estimation::PopulationEstimates::MLE
##   Estimation::IndividualEstimates::Estimates
##   Estimation::IndividualEstimates::RandomEffects
##   Estimation::Residuals::ResidualTable
##   Estimation::Predictions
##   Estimation::OFMeasures::Deviance
## 
## The following MESSAGEs were raised during the job execution:
##  estimation_successful: 1
##  covariance_step_run: 0
##  rounding_errors: 0
##  hessian_reset: 0
##  zero_gradients: 0
##  final_zero_gradients: 0
##  estimate_near_boundary: 0
##  s_matrix_singular: 0
##  significant_digits: 3
##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
## 
## Completed
## -- Mon Aug 15 13:58:41 2016

Once the model has run and the Standard Output (SO) object has returned it is possible to read that object in to R at any time.

# NM <-LoadSOObject("NONMEM/Simeoni_PAGE_Estimation_NM.SO.xml")

Results from NONMEM should be comparable to previous results

parameters_nm <- getPopulationParameters(NM,  what="estimates")$MLE
parameters_nm
##    POP_CL     POP_V    SD_ADD  OMEGA_CL   OMEGA_V 
## 0.8148350 8.0814500 0.4116460 0.0141337 0.0119764

ESTIMATE model parameters in WINBUGS

Assembling the new MOG for Winbugs. Note that we reuse the data and model. Two main changes are made: selecting the appropriate Task Properties Object for BUGS and exchanging the Prior Object instead of the Parameter Object.

BUGS.MOG <- createMogObj(dataObj = myDataObj, 
        priorObj = myPriorObj, 
        mdlObj = myModelObj, 
        taskObj = myTaskPropertiesObj_mcmc)

We can then write the MOG back out to an .mdl file.

mdlfile.BUGS <- paste0(model,"_BUGS.mdl")
writeMogObj(BUGS.MOG,mdlfile.BUGS)
#mdlfile.BUGS <- "Simeoni_PAGE_Estimation_BUGS.mdl"

There are two ways to run WinBUGS. The first is via the familiar estimate(…) statement. The “preprocessSteps = list(prepareWinbugsDs)” argument is mandatory.

# BUGS <- estimate(mdlfile, target="winbugs", subfolder="WinBUGS", preprocessSteps=list(prepareWinbugsDs))

Alternatively the runWinBUGS function is provided with the ddmore R package to allow the user some additional control over BUGS execution.

 BUGS<-runWinBUGS(mdlfile.BUGS,subfolder="WinBUGS")
## -- Mon Aug 15 13:58:43 2016
## New
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): rate column is
## missing. rate will be set to 0 by default.
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): ii column is
## missing. Default 0 values will be set.
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): cmt column is
## missing. cmt will be set to 1 by default.
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): cmt column is
## missing. cmt will be set to 1 by default.
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): addl column is
## missing. addl will be set to 0 by default.
## Warning in NMTRAN2BUGSdataconverter(model = model.pharmml): ss column is
## missing. ss will be set to 0 by default.
## Submitted
## Job 8fa28800-3ddf-40f7-a8ba-067b93ed6f2f progress:
## Running [ ...... ]
## Importing Results
## Copying the result data back to the local machine for job ID 8fa28800-3ddf-40f7-a8ba-067b93ed6f2f...
## From C:\Users\smith_mk\AppData\Local\Temp\4\RtmpYF4SFO\DDMORE.job35649cc6b71 to C:/SEE/MDL_IDE/workspace/InstallationTest/models/WinBUGS
## Done.
## 
## 
## The following main elements were parsed successfully:
##   RawResults
##   Estimation::PopulationEstimates::Bayesian
##   Estimation::PrecisionPopulationEstimates::Bayesian
## 
## The following MESSAGEs were raised during the job execution:
##  winbugs_message: success
## 
## Completed
## -- Mon Aug 15 14:00:47 2016
# BUGS <- LoadSOObject(file.path("WinBUGS","InstallationTest_DEQ_BUGS.SO.xml"))

getPopulationParameters(BUGS)$Bayesian
## Warning in getPopulationParameters(BUGS): Tried to fetch the parameter
## estimates, however section Estimation::IndividualEstimates not found in
## the SO Object.
##    Parameter        Mean   Median          SDP   Perc_2.5      Perc_5
## 1    lPOP_CL -0.20721360 -0.20770  0.035477621 -0.2750000 -0.26320000
## 2     lPOP_V  2.09322820  2.08900  0.037934270  2.0300000  2.03700000
## 3   OMEGA_CL  0.01312103  0.01160  0.006724577  0.0046820  0.00551105
## 4    OMEGA_V  0.01279900  0.01151  0.005940712  0.0053901  0.00604210
## 5     POP_CL  0.81335868  0.81240  0.028964354  0.7596000  0.76860000
## 6      POP_V  8.11685060  8.07900  0.309496561  7.6140000  7.66900000
## 7     SD_ADD  0.41600615  0.41140  0.043907612  0.3405000  0.35080000
## 8        TAU  5.96845610  5.90850  1.232344813  3.8300750  4.11305000
## 9     TAU_CL 95.20904200 86.18500 46.817657615 32.9710000 38.63000000
## 10     TAU_V 93.08502700 86.87500 38.585257038 35.4700000 41.88150000
##        Perc_25  Perc_50     Perc_75     Perc_95    Perc_97.5
## 1  -0.23060000 -0.20770  -0.1844000  -0.1486250  -0.13370250
## 2   2.06700000  2.08900   2.1180000   2.1610000   2.17100000
## 3   0.00856800  0.01160   0.0159875   0.0258895   0.03032925
## 4   0.00876925  0.01151   0.0152900   0.0238790   0.02819975
## 5   0.79410000  0.81240   0.8316000   0.8619750   0.87480000
## 6   7.90100000  8.07900   8.3170000   8.6800000   8.77100000
## 7   0.38520000  0.41140   0.4437000   0.4930950   0.51099500
## 8   5.07900000  5.90850   6.7410000   8.1270000   8.62600000
## 9  62.55250000 86.18500 116.7000000 181.4950000 213.60000000
## 10 65.41250000 86.87500 114.0000000 165.5000000 185.49750000
getPopulationParameters(BUGS, what="estimates")$Bayesian # o other options,  what= precisions, intervals, all
##          Parameter        Mean   Median
## lPOP_CL    lPOP_CL -0.20721360 -0.20770
## lPOP_V      lPOP_V  2.09322820  2.08900
## OMEGA_CL  OMEGA_CL  0.01312103  0.01160
## OMEGA_V    OMEGA_V  0.01279900  0.01151
## POP_CL      POP_CL  0.81335868  0.81240
## POP_V        POP_V  8.11685060  8.07900
## SD_ADD      SD_ADD  0.41600615  0.41140
## TAU            TAU  5.96845610  5.90850
## TAU_CL      TAU_CL 95.20904200 86.18500
## TAU_V        TAU_V 93.08502700 86.87500
# it is possible also to omitt $Bayesian
getIndividualParameters(BUGS)
## Warning in getIndividualParameters(BUGS): No values found for
## Estimates::Mean or RandomEffects::EffectMean in
## Estimation::IndividualEstimates slot of the SOObject.
## data frame with 0 columns and 0 rows

Here we use the coda package to examine and assess MCMC convergence.

library(coda)

By default the parameters monitored by BUGS are all parameters specified with random variable distributions in the PRIOR_VARIABLE_DEFINITION. In this example these are specified on the log scale. The following code transforms them back to the natural scale. Note that in the public release, MCMC nodes for monitoring will be specified via the Task Properties Object.

parameters_BUGS<- getPopulationParameters(BUGS, what="estimates")$Bayesian
estPars <- intersect(rownames(parameters_BUGS), names(parameters_nm))

coda MCMC trace and density plots

#to read only one chain
BUGSoutputPath <- file.path(getwd(),"WinBUGS")
coda_out <- read.coda(output.file=file.path(BUGSoutputPath,"output1.txt"), index.file=file.path(BUGSoutputPath,"outputIndex.txt"), quiet=T)
coda_pars <- coda_out[,estPars]
#to read more chains (3 in this case) use read.coda.interactive()
#coda_out <- read.coda.interactive()
# Enter in the order:
  #outputIndex.txt
  #output1.txt  
  #output2.txt  
  #output2.txt
## summary statistics
summary(coda_pars)
## 
## Iterations = 4001:14000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean       SD  Naive SE Time-series SE
## OMEGA_CL 0.01312 0.006725 6.725e-05      0.0001615
## OMEGA_V  0.01280 0.005941 5.941e-05      0.0001593
## POP_CL   0.81336 0.028964 2.896e-04      0.0019632
## POP_V    8.11685 0.309497 3.095e-03      0.0444528
## SD_ADD   0.41601 0.043908 4.391e-04      0.0009558
## 
## 2. Quantiles for each variable:
## 
##              2.5%      25%     50%     75%   97.5%
## OMEGA_CL 0.004682 0.008568 0.01160 0.01598 0.03030
## OMEGA_V  0.005394 0.008770 0.01151 0.01529 0.02819
## POP_CL   0.759600 0.794100 0.81240 0.83160 0.87480
## POP_V    7.614000 7.901000 8.07900 8.31700 8.77100
## SD_ADD   0.340500 0.385200 0.41140 0.44370 0.51080
## trace and density plots
plot(coda_pars, ask=F)

plot of chunk MCMC_Diagnostics plot of chunk MCMC_Diagnostics

## gelman-rubin diagnostics (if at least 2 chains)
#gelman.diag(coda_out)
#gelman.plot(coda_out)

Compare estimated parameters in the 3 tools

estPars <- intersect(rownames(parameters_BUGS), names(parameters_nm))
cbind(BUGS=parameters_BUGS[estPars,"Median"],
        NM=parameters_nm[estPars],
        MLX=parameters_mlx[estPars]) 
##             BUGS        NM     MLX
## OMEGA_CL 0.01160 0.0141337 0.01359
## OMEGA_V  0.01151 0.0119764 0.01232
## POP_CL   0.81240 0.8148350 0.81344
## POP_V    8.07900 8.0814500 8.07709
## SD_ADD   0.41140 0.4116460 0.41413

VPC of model

Update the (initial) values in the MDL Parameter Object Extract the structural and variability parameters from the SO object following NONMEM estimation.

structuralPar <- getPopulationParameters(MLX, what="estimates",block='structural')$MLE
variabilityPar <- getPopulationParameters(MLX, what="estimates",block='variability')$MLE

Update the parameter object using the ddmore “updateParObj” function. This function updates an R object of (S4) class “parObj”. The user chooses which block to update, what items within that block, and what to replace those items with. NOTE: that updateParObj can only update attributes which ALREADY EXIST in the MDL Parameter Object for that item. This ensures that valid MDL is preserved.

myParObjUpdated <- updateParObj(myParameterObj,block="STRUCTURAL",
        item=names(structuralPar),
        with=list(value=structuralPar))

myParObjUpdated <- updateParObj(myParObjUpdated,block="VARIABILITY",
        item=names(variabilityPar),
        with=list(value=variabilityPar))

Check that the appropriate initial values have been updated to the MLE values from the previous fit.

#print(myParObjUpdated@STRUCTURAL)
#print(myParObjUpdated@VARIABILITY)

Assembling the new MOG. Note that we reuse the data, model and tasks from the previous NONMEM run.

VPC.MOG <- createMogObj(dataObj = myDataObj, 
        parObj = myParObjUpdated, 
        mdlObj = myModelObj, 
        taskObj = myTaskPropertiesObj_focei)

We can then write the MOG back out to an .mdl file.

mdlfile.VPC <- paste0(model,"_VPC.mdl")
writeMogObj(VPC.MOG,mdlfile.VPC)

Similarly as above, ddmore “VPC.PsN” function can be used to run a VPC using PsN as target tool

vpcFiles <- VPC.PsN(mdlfile.VPC,samples=200, seed=12345,
        vpcOptions ="-n_simulation=10 -auto_bin=5, -min_points_in_bin=2",
        subfolder="VPC", plot=FALSE) 
## -- Mon Aug 15 14:00:58 2016
## New
## Submitted
## Job 69f21098-e3ab-4267-a511-85a3a6a0a734 progress:
## Running [ .... ]
## Importing Results
## Copying the result data back to the local machine for job ID 69f21098-e3ab-4267-a511-85a3a6a0a734...
## From C:\Users\smith_mk\AppData\Local\Temp\4\RtmpYF4SFO\DDMORE.job356467c26e1a to C:/SEE/MDL_IDE/workspace/InstallationTest/models/VPC
## Done.
## 
## 
## The following main elements were parsed successfully:
##   RawResults
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
##   SimulationSimulationBlock
## 
## The following MESSAGEs were raised during the job execution:
##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
## 
## Completed
## -- Mon Aug 15 14:02:21 2016

To replay the visualisation using information from the VPC SO file

plot(xpose.VPC(vpc.info= file.path("./VPC","vpc_results.csv"),
                vpctab= file.path("./VPC","vpctab"),
                main="Installation Test VPC"))

plot of chunk VPC_using_NONMEM_estimates

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)

library(mlxR)

mlxR contains a function pharmml2mlxtran(pharmMLFile) that performs the appropriate translation if needed.

#model.mlxtran <- file.path(getwd(),"Monolix","InstallationTest_MLX","InstallationTest_MLX_model.txt")

In future, we will be able to read Design Object information from MDL directly and create inputs to simulx by specifying as.simulx(mdlfile). For now though, we create the relevant inputs for simulx by hand.

p<-parameters_nm
print(p)
##    POP_CL     POP_V    SD_ADD  OMEGA_CL   OMEGA_V 
## 0.8148350 8.0814500 0.4116460 0.0141337 0.0119764

adm : dose administration times

adm <- list( time  = 0, 
              amount = 100,
              target='CENTRAL')

f : what to sample (“observe”) and at what times?

f <- list( name='CONC', time=seq(0,24,by=0.1))
y <- list( name='Y', time=c(0.001, 1,2,4,8,12,24))

g : define group sizes and assign treatments.

g <- list(size=12, level="longitudinal", treatment=adm)

Call simulx. simulx can take PharmML as an input - so the model does NOT need to be translated to MLXTRAN.

res <- simulx( model     = myPharmML,
               parameter = p,
               group = g,
               output    = list(f,y))

Plot simulx output

ggplot() + 
        geom_line(data=res$CONC, aes(x=time, y=CONC, colour=id))+
        geom_point(data=res$Y, aes(x=time, y=Y, colour=id))

plot of chunk Plot_of_simulated_data_using_NONMEM_estimates

Design evaluation using PFIM

Update the (initial) values in the MDL Parameter Object for PFIM evaluation Note that due to current limitations in PFIM, evaluation and optimisation of multiple doses (as in the Simeoni estimation model) is not straightforward to encode. For this demonstration, we will examine one single dose at time 8 (days) so that the tumour has time to grow before treatment, so that we can see the effect of treatment.

Although the multiple dose design is not supported by PFIM, the model stays the same and we can update the parameter values used for evaluating the study design using the results of estimation using NONMEM and use simulation results to provide the prediction of the total tumour weight (WTOT) at Day 8.

Assembling the new MOG.

PFIM.MOG <- createMogObj(designObj = myDesignObj, 
        parObj = myParObjUpdated, 
        mdlObj = myModelObj, 
        taskObj = myTaskPropertiesObj_Evaluate)

We can then write the MOG back out to an .mdl file.

mdlfile.PFIM <- "InstallationTest_DEQ_Design.mdl"
writeMogObj(PFIM.MOG,mdlfile.PFIM)
as.PharmML(mdlfile.PFIM)
## [1] "C:\\SEE\\MDL_IDE\\workspace\\InstallationTest\\models\\InstallationTest_DEQ_Design.xml"

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.PFIM)

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\\InstallationTest\\models\\InstallationTest_DEQ_Design.xml -o C:/SEE/MDL_IDE/workspace/InstallationTest/models/PFIM"
cat(readLines(file.path(getwd(),"PFIM","stdout.out")),sep="\n")
## PFIM 4.0 
##  
## Project:  Generated from MDL. MOG ID: outputMog
##  
## Date:  Mon Aug 15 14:02:37 2016
##  
## 
##  
## **************************** INPUT SUMMARY ********************************
##  
## Differential Equations form of the model:  
##  
## function (t, y, p) 
## {
##     CL <- p[1]
##     V <- p[2]
##     K <- ((CL)/(V))
##     CONC <- ((y[1])/(V))
##     yd1 <- ((-(K)) * (y[1]))
##     return(list(c(yd1), CONC))
## }
## 
##  
## Design:  
## Sample times for response: A 
##                           times subjects
## 1 c(1e-04, 1, 4, 8, 12, 24, 36)       12
## 
##  
## Initial Conditions at time 0 : 
##  
## 100 
## 
##  
## Random effect model: Trand =  2
##  
## Variance error model response A : ( 0.41413 + 0 *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
##  
## FIM saved in FIM.txt
##  
##  
## ******************* FISHER INFORMATION MATRIX ******************
##  
##             [,1]      [,2]         [,3]         [,4]       [,5]
## [1,] 151.7678399 0.0337299 0.000000e+00 0.000000e+00   0.000000
## [2,]   0.0337299 1.6478993 0.000000e+00 0.000000e+00   0.000000
## [3,]   0.0000000 0.0000000 4.201945e+02 2.046349e-03   5.922476
## [4,]   0.0000000 0.0000000 2.046349e-03 4.815815e+02   1.450527
## [5,]   0.0000000 0.0000000 5.922476e+00 1.450527e+00 699.782652
## 
##  
## ************************** EXPECTED STANDARD ERRORS ************************
##  
## ------------------------ Fixed Effects Parameters -------------------------
##  
##       Beta   StdError      RSE  
## CL 0.81344 0.08117291 9.978967 %
## V  8.07709 0.77899677 9.644523 %
## 
##  
## ------------------------- Variance of Inter-Subject Random Effects ----------------------
##  
##       omega2   StdError      RSE  
## CL 0.1165762 0.04878662 41.84957 %
## V  0.1109955 0.04556868 41.05453 %
## 
##  
## ------------------------ Standard deviation of residual error ------------------------ 
##  
##              Sigma   StdError      RSE  
## sig.interA 0.41413 0.03780469 9.128701 %
## 
##  
## ******************************* DETERMINANT ********************************
##  
## 35410920662
##  
## ******************************** CRITERION *********************************
##  
## 128.7738
##  
## 
##  
## ******************* EIGENVALUES OF THE FISHER INFORMATION MATRIX ******************
##  
##         FixedEffects VarianceComponents
## min       481.571849           1.647892
## max       699.917684         420.069087
## max/min     1.453402         254.913044
## 
##  
## ******************* CORRELATION MATRIX ******************
##  
##              [,1]         [,2]          [,3]          [,4]         [,5]
## [1,]  1.000000000 -0.002132848  0.000000e+00  0.000000e+00  0.000000000
## [2,] -0.002132848  1.000000000  0.000000e+00  0.000000e+00  0.000000000
## [3,]  0.000000000  0.000000000  1.000000e+00  2.274253e-05 -0.010921881
## [4,]  0.000000000  0.000000000  2.274253e-05  1.000000e+00 -0.002498771
## [5,]  0.000000000  0.000000000 -1.092188e-02 -2.498771e-03  1.000000000
## 
## 
##  
## Time difference of 0.1751168 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 311
## Warning: PopED Warning: No PopED-specific settings could be retrieved from
## modelling steps at line 277

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
## 2  0.0001 12.3805718     1     1        100           0
## 1  1.0000 11.1945717     1     1        100           0
## 6  4.0000  8.2755086     1     1        100           0
## 7  8.0000  5.5315177     1     1        100           0
## 3 12.0000  3.6973786     1     1        100           0
## 4 24.0000  1.1041873     1     1        100           0
## 5 36.0000  0.3297552     1     1        100           0

evaluate initial design

FIM <- evaluate.fim(poped.db) 
FIM
##            [,1]        [,2]       [,3]         [,4]         [,5]
## [1,] 58.7035778  -0.1165006 -0.1022716   25.7985718    8.9105434
## [2,] -0.1165006  97.4142590 -0.3435996  -10.0825113  -11.7880764
## [3,] -0.1022716  -0.3435996  1.2782903   -5.1636064   -7.9228357
## [4,] 25.7985718 -10.0825113 -5.1636064 1835.4578332    0.6646242
## [5,]  8.9105434 -11.7880764 -7.9228357    0.6646242 2985.6144113
det(FIM)
## [1] 38594088504
get_rse(FIM,poped.db)
##   bpop[1]   bpop[2]   bpop[3]    D[1,1]    D[2,2] 
##  31.62100  12.47047  11.11337 173.35002 149.88622