# DDMoRe: Demonstration of interoperability

Authors: Letizia Carrara, Celine Sarr, Mike K Smith Date: 18 July 2016

This example script is intended to demonstate interoperability using the DDMoRe Standalone Execution Environment (SEE).

The following steps are implemented in this workflow:

- Initialisation of the R console
- Exploratory graphical analysis
- Parameter estimation with Monolix and model evaluation in Xpose
- Parameter estimation in NONMEM using FOCE and model evaluation in Xpose
- Parameter estimation in Winbugs and MCMC diagnostics using coda
- Updating parameter estimates in the MDL Parameter Object using MLE values from NONMEM
- Performing a Visual Predictive Check (VPC) in PsN
- Using MLE values from NONMEM to simulate new observed values using the simulx function in the mlxR package.
- Updating parameter estimates with NONMEM and simulx output

# Initialisation

Clear workspace and set working directory

    rm(list=ls(all=F))

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-full\\MDL_IDE\\workspace"

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-full"

The “~” directory is the default MDL-IDE workspace directory

    projectDir <- file.path("~","NovartisDemo","models")
    exampleDir <- file.path(projectDir, "simeoni_pop")
    setwd(exampleDir)
    getwd()

    ## [1] "C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop"

Nominate the MDL file for analysis. The MDL file below describes the tumour growth of individual xenograft mice following innoculation with tumour cells.

    mdlfile <- "Simeoni_PAGE_Estimation_POP.mdl"
    model <- tools::file_path_sans_ext(mdlfile)

The ddmore R package has been written to allow users to read, write and work with MDL files and objects. getMDLObjects() returns all MDL objects from that file as an R list.

    myMDLObj <- getMDLObjects(mdlfile)
    length(myMDLObj)

    ## [1] 7

    names(myMDLObj)

    ## [1] "simeoni2004_data"         "simeoni2004_par"
    ## [3] "simeoni2004_mdl"          "simeoni2004_NONMEM_task"
    ## [5] "simeoni2004_Monolix_task" "simeoni2004_BUGS_task"
    ## [7] "simeoni2004_pri"

    #str(myMDLObj)

Next, 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.

    myDataObj <- getDataObjects(mdlfile)[[1]]
    datafile <- myDataObj@SOURCE$scrfile$file

    myParObj <- getParameterObjects(mdlfile)[[1]]
    myModObj <- getModelObjects(mdlfile)[[1]]

You can also refer to objects by their name.

    myTaskPropertiesObj_NM <- getTaskPropertiesObjects(mdlfile)[["simeoni2004_NONMEM_task"]]
    myTaskPropertiesObj_MLX <- getTaskPropertiesObjects(mdlfile)[["simeoni2004_Monolix_task"]]

# 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)
    head(myData)

    ##   ID TIME        DV AMT CMT MDV     GRP
    ## 1  1    0        NA   0   3   1 Treated
    ## 2  1    8 0.5607700   0   3   0 Treated
    ## 3  1   10 0.8402284   0   3   0 Treated
    ## 4  1   13 1.5304333   0   3   0 Treated
    ## 5  1   15 3.1849627   0   3   0 Treated
    ## 6  1   18 5.9940457   0   3   0 Treated

Extract only observation records

    myEDAData<-myData[myData$MDV==0,]

    library(ggplot2)

Plot the data using ggplot2 library

    plot1 <- ggplot() +
            geom_line(aes(y=DV,x=TIME,  color=as.factor(GRP), group=ID), data=myEDAData) +
            geom_point(aes(y=DV,x=TIME,  color=as.factor(GRP), group=ID), data=myEDAData) +
            labs(y=" Tumour volume (mm3)" ,x="Time (days)")
    print(plot1)

    plot2 <- ggplot(aes(y=DV,x=TIME, group=GRP), data=myEDAData) +
            geom_line()+
            geom_point() +
            facet_wrap(~ID) +
            labs(y=" Tumour volume (mm3)" ,x="Time (days)")
    print(plot2)

# 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 MDL file to do specific tasks.

To define the new MOG for estimation in Monolix using the ddmore createMogObj function:

    MLX.MOG <- createMogObj(dataObj = myDataObj,
            parObj = myParObj,
            mdlObj = myModObj,
            taskObj = myTaskPropertiesObj_MLX)

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

We can then write the MOG back out to an .mdl file using the ddmore writeMogObj function.

    writeMogObj(MLX.MOG,mdlfile.MLX)

The ddmore “estimate” function translates the contents of an .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")

    ## -- Tue Jul 19 13:54:26 2016
    ## New
    ## Submitted
    ## Job 1cabf21f-eb3b-4287-97b7-83992b52b424 progress:
    ## Running [ ....... ]
    ## Importing Results
    ## Copying the result data back to the local machine for job ID 1cabf21f-eb3b-4287-97b7-83992b52b424...
    ## From C:\Users\smith_mk\AppData\Local\Temp\1\RtmpgxsXuo\DDMORE.job319869eb5d96 to C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop/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
    ## -- Tue Jul 19 13:56:49 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/UseCase1.SO.xml")

## 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.

    parameters_mlx <- getPopulationParameters(mlx, what="estimates")$MLE
    parameters_mlx

    ##   LAMBDA0_POP   LAMBDA1_POP       PSI_POP       K10_POP       K21_POP
    ##       0.28103       0.75631      20.00000      20.83200       2.01100
    ##       K12_POP        K1_POP        K2_POP        W0_POP        V1_POP
    ##       0.14400       1.04245       0.67960       0.05199       0.81000
    ## OMEGA_LAMBDA0 OMEGA_LAMBDA1     omega_PSI     omega_K10     omega_K21
    ##       0.04044       0.35145       0.00000       0.00000       0.00000
    ##     omega_K12      omega_K1      omega_K2      omega_W0      omega_V1
    ##       0.00000       0.00000       0.00000       0.00000       0.00000
    ##             b
    ##       0.09579

Currently, we need a workaround to rename the “b” parameter from Monolix (b=proportional error parameter) to “CV” to match NONMEM.

    parameters_mlx <- reshape::rename(parameters_mlx, c(b="CV"))

We can also return the MLEs with associated precision / variances

    getPopulationParameters(mlx, what="precisions")

    ## $MLE
    ##        Parameter      MLE      SE   RSE
    ## 1              b  0.09579 0.00691  7.21
    ## 2         K1_POP  1.04245 0.00818  0.79
    ## 3        K10_POP 20.83200 0.00000  0.00
    ## 4        K12_POP  0.14400 0.00000  0.00
    ## 5         K2_POP  0.67960 0.00342  0.50
    ## 6        K21_POP  2.01100 0.00000  0.00
    ## 7    LAMBDA0_POP  0.28103 0.00319  1.13
    ## 8    LAMBDA1_POP  0.75631 0.06939  9.18
    ## 9       omega_K1  0.00000 0.00000  0.00
    ## 10     omega_K10  0.00000 0.00000  0.00
    ## 11     omega_K12  0.00000 0.00000  0.00
    ## 12      omega_K2  0.00000 0.00000  0.00
    ## 13     omega_K21  0.00000 0.00000  0.00
    ## 14 OMEGA_LAMBDA0  0.04044 0.00776 19.20
    ## 15 OMEGA_LAMBDA1  0.35145 0.06752 19.21
    ## 16     omega_PSI  0.00000 0.00000  0.00
    ## 17      omega_V1  0.00000 0.00000  0.00
    ## 18      omega_W0  0.00000 0.00000  0.00
    ## 19       PSI_POP 20.00000 0.00000  0.00
    ## 20        V1_POP  0.81000 0.00000  0.00
    ## 21        W0_POP  0.05199 0.00068  1.31

We can also return information about the model fit including LogLikelihood, AIC, BIC

    getEstimationInfo(mlx)

    ## $OFMeasures
    ## $OFMeasures$LogLikelihood
    ## $OFMeasures$LogLikelihood[[1]]
    ## [1] 7.645
    ##
    ##
    ## $OFMeasures$IndividualContribToLL
    ##    Subject ICtoLL
    ## 1        1  -7.17
    ## 2        2  -1.26
    ## 3        3  -5.55
    ## 4        4  -2.04
    ## 5        5  -4.21
    ## 6        6  -0.73
    ## 7        7  -3.41
    ## 8        8  -4.12
    ## 9        9   7.44
    ## 10      10   5.36
    ## 11      11   7.09
    ## 12      12   3.08
    ## 13      13   5.43
    ## 14      14   5.92
    ## 15      15  -3.04
    ## 16      16   4.86
    ##
    ## $OFMeasures$InformationCriteria
    ## $OFMeasures$InformationCriteria$AIC
    ## [1] 0.71
    ##
    ## $OFMeasures$InformationCriteria$BIC
    ## [1] 6.89
    ##
    ##
    ##
    ## $Messages
    ## list()

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.

    mlx.xpdb<-as.xpdb(mlx,datafile)

    ##
    ## Removed dose rows in rawData slot of SO to enable merge with Predictions data.

    basic.gof(mlx.xpdb)

    ind.plots(mlx.xpdb)

We can export graphs to PDF file

    pdf("GOF_MLX.pdf")
    print(basic.gof(mlx.xpdb))
    print(ind.plots(mlx.xpdb))
    dev.off()

## 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. Note that the NONMEM estimation method specified is FOCE. We will look at how to change this via R later.

    NM.MOG <- createMogObj(dataObj = myDataObj,
            parObj = myParObj,
            mdlObj = myModObj,
            taskObj = myTaskPropertiesObj_NM)

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

    ## -- Tue Jul 19 13:56:52 2016
    ## New
    ## Submitted
    ## Job 33260fa5-45d7-479a-80df-841d9d5b430b progress:
    ## Running [ ....... ]
    ## Importing Results
    ## Copying the result data back to the local machine for job ID 33260fa5-45d7-479a-80df-841d9d5b430b...
    ## From C:\Users\smith_mk\AppData\Local\Temp\1\RtmpgxsXuo\DDMORE.job319860cd426e to C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop/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.8
    ##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
    ##
    ## Completed
    ## -- Tue Jul 19 13:59:15 2016

Load previous results

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

Results from NONMEM should be comparable with results from MONOLIX

    parameters_nm <- getPopulationParameters(NM, what="estimates")$MLE

We can compare the results across tools. We sort on the names of parameters from NONMEM output

    estPars <- names(parameters_nm)
    cbind(  NM=parameters_nm[estPars],
            MLX=parameters_mlx[estPars])

    ##                        NM      MLX
    ## LAMBDA0_POP    0.26595100  0.28103
    ## LAMBDA1_POP    0.77905900  0.75631
    ## CV             0.09886030  0.09579
    ## K1_POP         1.07391000  1.04245
    ## K2_POP         0.61950600  0.67960
    ## W0_POP         0.05950020  0.05199
    ## K10_POP       20.83200000 20.83200
    ## K12_POP        0.14400000  0.14400
    ## K21_POP        2.01100000  2.01100
    ## V1_POP         0.81000000  0.81000
    ## PSI_POP       20.00000000 20.00000
    ## OMEGA_LAMBDA0  0.00113818  0.04044
    ## OMEGA_LAMBDA1  0.13398600  0.35145

Covariance step cannot be requested in the current version of the Framework

    #getPopulationParameters(NM, what="precisions")

As before we can request additional information about the NONMEM estimation process.

    print(getEstimationInfo(NM))

    ## $OFMeasures
    ## $OFMeasures$Deviance
    ## $OFMeasures$Deviance[[1]]
    ## [1] -254.3298
    ##
    ##
    ##
    ## $Messages
    ## $Messages$Info
    ## $Messages$Info$estimation_successful
    ## [1] "1"
    ##
    ## $Messages$Info$covariance_step_run
    ## [1] "0"
    ##
    ## $Messages$Info$rounding_errors
    ## [1] "0"
    ##
    ## $Messages$Info$hessian_reset
    ## [1] "0"
    ##
    ## $Messages$Info$zero_gradients
    ## [1] "0"
    ##
    ## $Messages$Info$final_zero_gradients
    ## [1] "0"
    ##
    ## $Messages$Info$estimate_near_boundary
    ## [1] "0"
    ##
    ## $Messages$Info$s_matrix_singular
    ## [1] "0"
    ##
    ## $Messages$Info$significant_digits
    ## [1] "3.8"
    ##
    ## $Messages$Info$nmoutput2so_version
    ## [1] "This SOBlock was created with nmoutput2so version 4.5.27"

Again we convert the SO object to an Xpose format database and perform model diagnostics.

    nm.xpdb<-as.xpdb(NM,datafile)

    ##
    ## Removed dose rows in rawData+Predictions slot of SO to enable merge with Residuals data.
    ##
    ## Residuals data does not currently contain dose rows in output from Nonmem executions.

Perform some basic goofness of fit

    basic.gof(nm.xpdb)

    ind.plots(nm.xpdb)

Export graphs to PDF file

    pdf("GOF_NM.pdf")
    print(basic.gof(mlx.xpdb))
    print(ind.plots(mlx.xpdb))
    dev.off()

## Change estimation algorithm to SAEM

Let's change the estimation algorithm from FOCE to SAEM.

    myTaskPropertiesObj_NM_SAEM <- myTaskPropertiesObj_NM
    myTaskPropertiesObj_NM_SAEM@ESTIMATE$algo <- "saem"

create new MOG, write to file and estimate.

    NM.MOG <- createMogObj(dataObj = myDataObj,
            parObj = myParObj,
            mdlObj = myModObj,
            taskObj = myTaskPropertiesObj_NM_SAEM)

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

    NM.SAEM <- estimate(mdlfile.NM.SAEM, target="NONMEM", subfolder="Nonmem_SAEM")

    ## -- Tue Jul 19 13:59:18 2016
    ## New
    ## Submitted
    ## Job 03345865-87cd-401f-95f9-1a4e7e04704f progress:
    ## Running [ ........................................................... ]
    ## Importing Results
    ## Copying the result data back to the local machine for job ID 03345865-87cd-401f-95f9-1a4e7e04704f...
    ## From C:\Users\smith_mk\AppData\Local\Temp\1\RtmpgxsXuo\DDMORE.job31983ecf61d6 to C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop/Nonmem_SAEM
    ## 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
    ##  estimate_near_boundary: 0
    ##  s_matrix_singular: 0
    ##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
    ##
    ## Completed
    ## -- Tue Jul 19 14:19:03 2016

    parameters_nm_saem <- getPopulationParameters(NM.SAEM, what="estimates")$MLE

We can compare the results across tools. We sort on the names of parameters from NONMEM output

    estPars <- names(parameters_nm)
    cbind(  NM=parameters_nm[estPars],
            NM.SAEM=parameters_nm_saem[estPars],
            MLX=parameters_mlx[estPars])

    ##                        NM     NM.SAEM      MLX
    ## LAMBDA0_POP    0.26595100 3.36109e-01  0.28103
    ## LAMBDA1_POP    0.77905900 6.44414e-01  0.75631
    ## CV             0.09886030 1.79577e-01  0.09579
    ## K1_POP         1.07391000 7.87623e-01  1.04245
    ## K2_POP         0.61950600 7.71906e-01  0.67960
    ## W0_POP         0.05950020 3.01733e-02  0.05199
    ## K10_POP       20.83200000 2.08320e+01 20.83200
    ## K12_POP        0.14400000 1.44000e-01  0.14400
    ## K21_POP        2.01100000 2.01100e+00  2.01100
    ## V1_POP         0.81000000 8.10000e-01  0.81000
    ## PSI_POP       20.00000000 2.00000e+01 20.00000
    ## OMEGA_LAMBDA0  0.00113818 1.77595e-09  0.04044
    ## OMEGA_LAMBDA1  0.13398600 1.24513e-01  0.35145

## 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(NM, what="estimates",block='structural')$MLE
    variabilityPar <- getPopulationParameters(NM, 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.

    myParObj <- getParameterObjects(mdlfile)[[1]]
    myParObjUpdated <- updateParObj(myParObj,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_ALL <- createMogObj(dataObj = myDataObj,
            parObj = myParObjUpdated,
            mdlObj = myModObj,
            taskObj = myTaskPropertiesObj_NM)

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

    mdlfile.VPC_ALL <- paste0(model,"_VPC_ALL.mdl")
    writeMogObj(VPC.MOG_ALL,mdlfile.VPC_ALL)

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

    vpcFiles <- VPC.PsN(mdlfile.VPC_ALL,samples=20, seed=12345,
            vpcOptions ="-n_simulation=10  -auto_bin=3  -stratify_on=GRP" ,
            subfolder="VPC", plot=FALSE)

    ## -- Tue Jul 19 14:19:09 2016
    ## New
    ## Submitted
    ## Job f37be0aa-4b00-4008-8a9f-e9cff82e4308 progress:
    ## Running [ ... ]
    ## Importing Results
    ## Copying the result data back to the local machine for job ID f37be0aa-4b00-4008-8a9f-e9cff82e4308...
    ## From C:\Users\smith_mk\AppData\Local\Temp\1\RtmpgxsXuo\DDMORE.job31985fd1111a to C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop/VPC
    ## Done.
    ##
    ##
    ## The following main elements were parsed successfully:
    ##   RawResults
    ##   SimulationSimulationBlock
    ##   SimulationSimulationBlock
    ##
    ## The following MESSAGEs were raised during the job execution:
    ##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
    ##
    ## Completed
    ## -- Tue Jul 19 14:20:12 2016

To replay the visualisation using information from the VPC SO file

    print(xpose.VPC(vpc.info= file.path("./VPC",vpcFiles@RawResults@DataFiles$PsN_VPC_results$path),
                    vpctab= file.path("./VPC",vpcFiles@RawResults@DataFiles$PsN_VPC_vpctab$path),
                    by="GRP",
                    main="VPC SIMEONI model"))

## 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)
    model.mlxtran <- pharmml2mlxtran(myPharmML, "InstallationTest_Model.txt")

    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. Use parameter values from the FOCE estimation

    p <- parameters_nm

adm : dose administration times

    adm <- list( time  = c(8, 12, 16),
            amount = 30,
            target = 'Q1')

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

    f <- list( name='WTOT', time=seq(0,45,by=0.1))
    y <- list( name='Y', time=c(0.001, 8, 12, 16, 20, 24, 28, 32, 35))

g : define group sizes and assign treatments.

    g0 <- list(size=8, level='longitudinal')
    g1 <- list(size=8, 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 = list(g0,g1),
            output    = list(f,y))

Plot simulx output

    ggplot() +
            geom_line(data=res$WTOT, aes(x=time, y=WTOT, colour=id))+
            geom_point(data=res$Y, aes(x=time, y=Y, colour=id))+
            scale_colour_discrete(name="Experimental\n condition",
                    breaks=c("1", "2"),
                    labels=c("Control", "Treatment"))

## ESTIMATE model parameters in WINBUGS

BUGS dataset conversion currently has problems with categorical covariates. So we read the Data Object in MDL and change the use of GRP variable to “ignore”.

    BUGSDataObj <- myDataObj
    BUGSDataObj@DATA_INPUT_VARIABLES$GRP$use <- "ignore"

    BUGSModObj <- myModObj
    slot(BUGSModObj,"COVARIATES") <- list()

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 = BUGSDataObj,
            priorObj = getPriorObjects(mdlfile)[[1]],
            mdlObj = BUGSModObj,
            taskObj = getTaskPropertiesObjects(mdlfile)[["simeoni2004_BUGS_task"]])

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"
    file.exists(mdlfile.BUGS)

    ## [1] TRUE

Winbugs converter and function call has not been integrated into the SEE yet, so we cannot run:

    # estimate(mdlfile.BUGS, target="BUGS")

in the meantime, an R script has been prepared to call the converter, generate Winbugs input. This involves conversion of the NONMEM format dataset to BUGS (R list) format and preparation of the model file, data file, initial values for the MCMC and creation of the script file for running BUGS in batch mode.

    source(file.path(ddmore:::DDMORE.checkConfiguration(),"BlackBoxWinBUGS","ddmorewinbugs.R"))

    BUGS<-run.winbugs(mdlfile.BUGS,subfolder="WinBUGS")

    ## Warning in normalizePath(path.expand(path), winslash, mustWork):
    ## path[1]="Simeoni_PAGE_Estimation_POP_BUGS.mdl": The system cannot find the
    ## file specified

    ## Error in ddmore:::execute(mdlfile, target = "winbugs", importSO = FALSE, : Illegal Argument: file  Simeoni_PAGE_Estimation_POP_BUGS.mdl  does not exist.

    getPopulationParameters(BUGS)$Bayesian

    ## Error in is(SOObject, class2 = "StandardOutputObject"): object 'BUGS' not found

    getPopulationParameters(BUGS, what="estimates")$Bayesian # o other options,  what= precisions, intervals, all

    ## Error in is(SOObject, class2 = "StandardOutputObject"): object 'BUGS' not found

    parameters_BUGS<- getPopulationParameters(BUGS, what="estimates")$Bayesian

    ## Error in is(SOObject, class2 = "StandardOutputObject"): object 'BUGS' not found

    # it is possible also to omitt $Bayesian
    getIndividualParameters(BUGS)

    ## Error in getIndividualParameters(BUGS): object 'BUGS' not found

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

    library(coda)

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

    ## Warning in file(file, "rt"): cannot open file
    ## 'C:/see-full/MDL_IDE/workspace/NovartisDemo/models/simeoni_pop/WinBUGS/outputIndex.txt':
    ## No such file or directory

    ## Error in file(file, "rt"): cannot open the connection

Summary statistics of MCMC samples

    summary(coda_out)

    ## Error in summary(coda_out): object 'coda_out' not found

    ## trace and density plots
    #pdf("MCMC_diagnostics_BUGS.pdf")
    #plot(coda_out, ask=F)
    #dev.off()
    ## gelman-rubin diagnostics (if at least 2 chains)
    #gelman.diag(coda_out)
    #gelman.plot(coda_out)

The SO conversion for BUGS output is not yet complete, but eventually we should be able to perform Xpose diagnostics after estimation with BUGS as we have with Monolix and NONMEM.

    # bugs<-as.xpdb(BUGS,datafile)

## Compare estimated parameters in the 3 tools

Note that in this example, we may need to run the MCMC sampler in BUGS to get better estimates of model parameters (see the MCMC diagnostics above), but for expediency we have limited the sampler to 5000 samples from the posterior distribution.

    estPars <- intersect(rownames(parameters_BUGS), names(parameters_nm))

    ## Error in rownames(parameters_BUGS): object 'parameters_BUGS' not found

    cbind(BUGS=parameters_BUGS[estPars,"Median"],
            NM=parameters_nm[estPars],
            MLX=parameters_mlx[estPars])

    ## Error in cbind(BUGS = parameters_BUGS[estPars, "Median"], NM = parameters_nm[estPars], : object 'parameters_BUGS' not found