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:
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"]]
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)
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"
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.
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))
mlx.xpdb@Data$ID <- as.numeric(as.character(mlx.xpdb@Data$ID))
print(ind.plots(mlx.xpdb))
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
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)
## gelman-rubin diagnostics (if at least 2 chains)
#gelman.diag(coda_out)
#gelman.plot(coda_out)
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
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"))
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))
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
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)
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