UseCase1_1 : ODE model for warfarin population pharmacokinetics (extended task properties)

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 wiht the coursor 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=F))
mydir <- file.path(Sys.getenv("MDLIDE_WORKSPACE_HOME"),"UseCasesDemo")
setwd(mydir)

Set name of .mdl file and dataset for future tasks

uc <- "UseCase1_1"
datafile <- "warfarin_conc.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))
cache.path <- file.path(mydir,uc, 'cache')
dir.create(cache.path)
opts_chunk$set(cache.path=cache.path)
figure.path <- file.path(mydir,uc, 'figure')
dir.create(figure.path)
opts_chunk$set(figure.path=figure.path)

List files available in working directory

list.files()
## [1] "cache"             "figure"            "UseCase1_1.mdl"   
## [4] "warfarin_conc.csv"

Introduction to 'ddmore' R package

View objects within the .mdl file

Use 'ddmore' function getMDLObjects() to retrieve model object(s) from an existing .mdl file. This function reads the MDL in an .mdl file and parses the MDL code for each MDL Object into an R list of objects of appropriate types with names corresponding to the MDL Object names given in the file.

myMDLObj <- getMDLObjects(mdlfile)
length(myMDLObj)
## [1] 6
names(myMDLObj)
## [1] "warfarin_PK_ODE_dat" "warfarin_PK_ODE_par" "warfarin_PK_ODE_mdl"
## [4] "SAEM_task"           "BUGS_task"           "FOCEI_task"

Use 'ddmore' function getDataObjects() to retrieve only data object(s) from an existing .mdl file. This function returns a list of Parameter Object(s) from which we select the first element Hover over the variable name to see its structure

myDataObj <- getDataObjects(mdlfile)[[1]]

Use 'ddmore' function getParameterObjects() to retrieve only parameter object(s) from an existing .mdl file

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

Use 'ddmore' function getModelObjects() to retrieve only model object(s) from an existing .mdl file.

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

Use 'ddmore' function getTaskPropertiesObjects() to retrieve only task properties object(s) from an existing .mdl file

myTaskObj <- getTaskPropertiesObjects(mdlfile)[[1]]

Objects can also be retrieved by name

myTaskObj.FOCEI <- getTaskPropertiesObjects(mdlfile)[["FOCEI_task"]]
myTaskObj.SAEM <- getTaskPropertiesObjects(mdlfile)[["SAEM_task"]]

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 16:36:34 2016
## New
## Submitted
## Job c2fd6c99-b3c5-4971-9944-d99010a8d21a progress:
## Running [ .......... ]
## Importing Results
## Copying the result data back to the local machine for job ID c2fd6c99-b3c5-4971-9944-d99010a8d21a...
## From C:\Users\zparra\AppData\Local\Temp\Rtmpk5tG0t\DDMORE.job1c842dea41b7 to D:/SEE-Prod5_RC4/MDL_IDE/workspace/UseCasesDemo/UseCase1_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::Residuals::ResidualTable
##   Estimation::Predictions
##   Estimation::OFMeasures::IndividualContribToLL
##   Estimation::OFMeasures::InformationCriteria
##   Estimation::OFMeasures::LogLikelihood
## 
## Completed
## -- Tue Aug 16 16:40:01 2016

Checking that the Monolix task settings have been provided to the MLXTRAN Project file correctly.

myTaskObj.SAEM
## An object of class "taskObj"
## Slot "ESTIMATE":
## $algo
## [1] "saem"
## 
## $blocks
## $blocks[[1]]
## $blocks[[1]]$TARGET_SETTINGS
## $blocks[[1]]$TARGET_SETTINGS$INTER
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$NBURN
## [1] "3000"
## 
## $blocks[[1]]$TARGET_SETTINGS$NITER
## [1] "500"
## 
## $blocks[[1]]$TARGET_SETTINGS$PRINT
## [1] "100"
## 
## $blocks[[1]]$TARGET_SETTINGS$SEED
## [1] "1556678"
## 
## $blocks[[1]]$TARGET_SETTINGS$ISAMPLE
## [1] "2"
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_MATRIX
## [1] "\"R\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_PRINT
## [1] "\"E\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_UNCONDITIONAL
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_SIGL
## [1] "12"
## 
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs$target
## [1] "\"NONMEM\""
## 
## 
## 
## $blocks[[1]]$.subtype
## [1] "BlockStmt"
## 
## 
## $blocks[[2]]
## $blocks[[2]]$TARGET_SETTINGS
## $blocks[[2]]$TARGET_SETTINGS$graphicsSettings
## [1] "\"tables.xmlx\""
## 
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs$target
## [1] "\"MONOLIX\""
## 
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs$settingsFile
## [1] "[\"tables.xmlx\"]"
## 
## 
## 
## $blocks[[2]]$.subtype
## [1] "BlockStmt"
## 
## 
## 
## 
## Slot "SIMULATE":
## NULL
## 
## Slot "EVALUATE":
## NULL
## 
## Slot "OPTIMISE":
## NULL
## 
## Slot "name":
## [1] "SAEM_task"
scan(file.path(getwd(),"Monolix",uc,paste(uc,"_project.mlxtran",sep="")), strip.white=T, what="character", sep="\n")
##  [1] "DATA:"                                                            
##  [2] "path = \"%MLXPROJECT%/\","                                        
##  [3] "file  =\"warfarin_conc.csv\","                                    
##  [4] "headers = {ID,TIME,IGNORE,DOSE,YTYPE,Y,MDV,COV}"                  
##  [5] "columnDelimiter = \",\""                                          
##  [6] "VARIABLES:"                                                       
##  [7] "logtWT [use=cov],"                                                
##  [8] "INDIVIDUAL:"                                                      
##  [9] "CL = { distribution=lognormal, covariate={logtWT} , iiv=yes },"   
## [10] "V = { distribution=lognormal, covariate={logtWT} , iiv=yes },"    
## [11] "KA = { distribution=lognormal, iiv=yes },"                        
## [12] "TLAG = { distribution=lognormal, iiv=yes },"                      
## [13] "CORRELATION:"                                                     
## [14] "correlationIIV={{CL,V}}"                                          
## [15] "OBSERVATION:"                                                     
## [16] "Y = { type=continuous, prediction=CC, error=combined1},"          
## [17] "STRUCTURAL_MODEL:"                                                
## [18] "file=\"mlxt:UseCase1_1_model\","                                  
## [19] "path=\"%MLXPROJECT%\","                                           
## [20] "output={CC}"                                                      
## [21] "TASKS:"                                                           
## [22] "globalSettings={"                                                 
## [23] "withVariance=no,"                                                 
## [24] "settingsGraphics=\"%MLXPROJECT%/tables.xmlx\","                   
## [25] "},"                                                               
## [26] "; workflow"                                                       
## [27] "estimatePopulationParameters("                                    
## [28] "initialValues={"                                                  
## [29] "beta_{CL,logtWT} = 0.75 [method=FIXED],"                          
## [30] "beta_{V,logtWT} = 1 [method=FIXED],"                              
## [31] "omega_{CL} = 0.1,"                                                
## [32] "omega_{KA} = 0.1,"                                                
## [33] "omega_{TLAG} = 0.1 [method=FIXED],"                               
## [34] "omega_{V} = 0.1,"                                                 
## [35] "pop_{CL} = 0.1,"                                                  
## [36] "pop_{KA} = 0.362,"                                                
## [37] "pop_{TLAG} = 1,"                                                  
## [38] "pop_{V} = 8,"                                                     
## [39] "a_{Y} = 0.1,"                                                     
## [40] "b_{Y} = 0.1,"                                                     
## [41] "})"                                                               
## [42] "estimateFisherInformationMatrix( method={linearization} ),"       
## [43] "estimateIndividualParameters( method={conditionalDistribution} ),"
## [44] "estimateLogLikelihood(method={linearization}),"                   
## [45] "displayGraphics(),"

SAEM Estimation with NONMEM

NM <- estimate(mdlfile, target="NONMEM", subfolder="NONMEM")
## -- Tue Aug 16 16:40:01 2016
## New
## Submitted
## Job 4cdc34d2-2e40-4027-97b7-27fc267ee841 progress:
## Running [ ............................................................................................................. ]
## Importing Results
## Copying the result data back to the local machine for job ID 4cdc34d2-2e40-4027-97b7-27fc267ee841...
## From C:\Users\zparra\AppData\Local\Temp\Rtmpk5tG0t\DDMORE.job1c843cad369d to D:/SEE-Prod5_RC4/MDL_IDE/workspace/UseCasesDemo/UseCase1_1/NONMEM
## Done.
## 
## 
## The following main elements were parsed successfully:
##   RawResults
##   Estimation::PopulationEstimates::MLE
##   Estimation::PrecisionPopulationEstimates::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: 1
##  covariance_step_successful: 1
##  covariance_step_warnings: 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 Aug 16 17:16:39 2016
myTaskObj.SAEM
## An object of class "taskObj"
## Slot "ESTIMATE":
## $algo
## [1] "saem"
## 
## $blocks
## $blocks[[1]]
## $blocks[[1]]$TARGET_SETTINGS
## $blocks[[1]]$TARGET_SETTINGS$INTER
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$NBURN
## [1] "3000"
## 
## $blocks[[1]]$TARGET_SETTINGS$NITER
## [1] "500"
## 
## $blocks[[1]]$TARGET_SETTINGS$PRINT
## [1] "100"
## 
## $blocks[[1]]$TARGET_SETTINGS$SEED
## [1] "1556678"
## 
## $blocks[[1]]$TARGET_SETTINGS$ISAMPLE
## [1] "2"
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_MATRIX
## [1] "\"R\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_PRINT
## [1] "\"E\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_UNCONDITIONAL
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_SIGL
## [1] "12"
## 
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs$target
## [1] "\"NONMEM\""
## 
## 
## 
## $blocks[[1]]$.subtype
## [1] "BlockStmt"
## 
## 
## $blocks[[2]]
## $blocks[[2]]$TARGET_SETTINGS
## $blocks[[2]]$TARGET_SETTINGS$graphicsSettings
## [1] "\"tables.xmlx\""
## 
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs$target
## [1] "\"MONOLIX\""
## 
## $blocks[[2]]$TARGET_SETTINGS$blkAttrs$settingsFile
## [1] "[\"tables.xmlx\"]"
## 
## 
## 
## $blocks[[2]]$.subtype
## [1] "BlockStmt"
## 
## 
## 
## 
## Slot "SIMULATE":
## NULL
## 
## Slot "EVALUATE":
## NULL
## 
## Slot "OPTIMISE":
## NULL
## 
## Slot "name":
## [1] "SAEM_task"
scan(file.path(getwd(),"NONMEM",paste(uc,".ctl",sep="")), what="character", sep="\n", strip.white=T)
##  [1] "; Script generated by the pharmML2Nmtran Converter v.0.4.0"                                         
##  [2] "; Source \t: PharmML 0.8.1"                                                                          
##  [3] "; Target \t: NMTRAN 7.3.0"                                                                           
##  [4] "; Model \t: UseCase1_1"                                                                              
##  [5] "; Model Name: 10May2016 Task Properties check"                                                      
##  [6] "; Dated \t: Tue Aug 16 16:40:12 CEST 2016"                                                           
##  [7] "$PROBLEM  my Problem Statement"                                                                     
##  [8] "$INPUT  ID TIME WT=DROP AMT DVID DV MDV LOGTWT"                                                     
##  [9] "$DATA \"warfarin_conc.csv\" IGNORE=@"                                                               
## [10] "$SUBS ADVAN13 TOL=6"                                                                                
## [11] "$MODEL"                                                                                             
## [12] "COMP (COMP1 DEFDOSE) \t;GUT"                                                                         
## [13] "COMP (COMP2) \t;CENTRAL"                                                                             
## [14] "$PK"                                                                                                
## [15] "POP_CL = THETA(1)"                                                                                  
## [16] "POP_V = THETA(2)"                                                                                   
## [17] "POP_KA = THETA(3)"                                                                                  
## [18] "POP_TLAG = THETA(4)"                                                                                
## [19] "RUV_PROP = THETA(5)"                                                                                
## [20] "RUV_ADD = THETA(6)"                                                                                 
## [21] "BETA_CL_WT = THETA(7)"                                                                              
## [22] "BETA_V_WT = THETA(8)"                                                                               
## [23] "ETA_CL = ETA(1)"                                                                                    
## [24] "ETA_V = ETA(2)"                                                                                     
## [25] "ETA_KA = ETA(3)"                                                                                    
## [26] "ETA_TLAG = ETA(4)"                                                                                  
## [27] "PPV_TLAG = 0.0"                                                                                     
## [28] "MU_1 = LOG(POP_CL) + BETA_CL_WT * LOGTWT"                                                           
## [29] "CL =  EXP(MU_1 + ETA(1)) ;"                                                                         
## [30] "MU_2 = LOG(POP_V) + BETA_V_WT * LOGTWT"                                                             
## [31] "V =  EXP(MU_2 + ETA(2)) ;"                                                                          
## [32] "MU_3 = LOG(POP_KA)"                                                                                 
## [33] "KA =  EXP(MU_3 + ETA(3)) ;"                                                                         
## [34] "MU_4 = LOG(POP_TLAG)"                                                                               
## [35] "TLAG =  EXP(MU_4 + ETA(4)) ;"                                                                       
## [36] "A_0(1) = 0"                                                                                         
## [37] "A_0(2) = 0"                                                                                         
## [38] "$DES"                                                                                               
## [39] "GUT_DES = A(1)"                                                                                     
## [40] "CENTRAL_DES = A(2)"                                                                                 
## [41] "RATEIN_DES = 0"                                                                                     
## [42] "IF (T.GE.TLAG) THEN"                                                                                
## [43] "RATEIN_DES = (GUT_DES*KA)"                                                                          
## [44] "ENDIF"                                                                                              
## [45] "CC_DES = (CENTRAL_DES/V)"                                                                           
## [46] "DADT(1) = -(RATEIN_DES)"                                                                            
## [47] "DADT(2) = (RATEIN_DES-((CL*CENTRAL_DES)/V))"                                                        
## [48] "$ERROR"                                                                                             
## [49] "EPS_Y = EPS(1)"                                                                                     
## [50] "GUT = A(1)"                                                                                         
## [51] "CENTRAL = A(2)"                                                                                     
## [52] "RATEIN = 0"                                                                                         
## [53] "IF (TIME.GE.TLAG) THEN"                                                                             
## [54] "RATEIN = (GUT*KA)"                                                                                  
## [55] "ENDIF"                                                                                              
## [56] "CC = (CENTRAL/V)"                                                                                   
## [57] "IPRED = CC"                                                                                         
## [58] "W = RUV_ADD+RUV_PROP*IPRED"                                                                         
## [59] "Y = IPRED+W*EPS_Y"                                                                                  
## [60] "IRES = DV - IPRED"                                                                                  
## [61] "IWRES = IRES/W"                                                                                     
## [62] "$THETA"                                                                                             
## [63] "( 0.001 , 0.1 )\t;POP_CL"                                                                            
## [64] "( 0.001 , 8.0 )\t;POP_V"                                                                             
## [65] "( 0.001 , 0.362 )\t;POP_KA"                                                                          
## [66] "( 0.001 , 1.0 )\t;POP_TLAG"                                                                          
## [67] "( 0.0 , 0.1 )\t;RUV_PROP"                                                                            
## [68] "( 1.0E-4 , 0.1 )\t;RUV_ADD"                                                                          
## [69] "(0.75  FIX )\t;BETA_CL_WT"                                                                           
## [70] "(1.0  FIX )\t;BETA_V_WT"                                                                             
## [71] "$OMEGA BLOCK(2) CORRELATION SD"                                                                     
## [72] "(0.1 )\t;PPV_CL"                                                                                     
## [73] "(0.01 )\t;CORR_CL_V"                                                                                 
## [74] "(0.1 )\t;PPV_V"                                                                                      
## [75] "$OMEGA"                                                                                             
## [76] "(0.1 SD )\t;PPV_KA"                                                                                  
## [77] "(0.1  FIX SD )\t;PPV_TLAG"                                                                           
## [78] "$SIGMA"                                                                                             
## [79] "1.0 FIX"                                                                                            
## [80] "$EST METHOD=SAEM  INTER NBURN=3000 NITER=500 PRINT=100 SEED=1556678 ISAMPLE=2"                      
## [81] "$COV  MATRIX=R PRINT=E UNCONDITIONAL SIGL=12"                                                       
## [82] "$TABLE  ID TIME AMT DVID MDV LOGTWT PRED RES WRES DV IPRED IRES IWRES Y NOAPPEND NOPRINT FILE=sdtab"
## [83] "$TABLE  ID CL V KA TLAG ETA_CL ETA_V ETA_KA ETA_TLAG NOAPPEND NOPRINT FILE=patab"                   
## [84] "$TABLE  ID LOGTWT NOAPPEND NOPRINT FILE=cotab"

Change estimation method to FOCEI (for speed)

Assembling the new Modelling Object Group (MOG). Note that we reuse the data, parameters and model from the MOG and only task properties is changed

myNewerMOG <- createMogObj(dataObj = getDataObjects(mdlfile)[[1]], 
        parObj = getParameterObjects(mdlfile)[[1]], 
        mdlObj = getModelObjects(mdlfile)[[1]], 
        taskObj = getTaskPropertiesObjects(mdlfile)[["FOCEI_task"]],
        info=list(problemStmt="Warfarin Pop PK base model", name="UseCase1 variant 1"))

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

mdlfile.FOCEI <- paste0(uc,"_FOCEI.mdl")
writeMogObj(myNewerMOG,mdlfile.FOCEI)

Test estimation using this new MOG in NONMEM via PsN

NM.FOCEI <- estimate(mdlfile.FOCEI, target="PsN", subfolder="NONMEM_FOCEI")
## -- Tue Aug 16 17:16:48 2016
## New
## Submitted
## Job f2e8348e-21da-4ef1-bca6-408ce0c53123 progress:
## Running [ ............... ]
## Importing Results
## Copying the result data back to the local machine for job ID f2e8348e-21da-4ef1-bca6-408ce0c53123...
## From C:\Users\zparra\AppData\Local\Temp\Rtmpk5tG0t\DDMORE.job1c8433f5390a to D:/SEE-Prod5_RC4/MDL_IDE/workspace/UseCasesDemo/UseCase1_1/NONMEM_FOCEI
## 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 WARNINGs were raised during the job execution:
##  covariance_step_successful: 0
##  covariance_step_warnings: 1
## 
## The following MESSAGEs were raised during the job execution:
##  estimation_successful: 1
##  covariance_step_run: 1
##  rounding_errors: 0
##  hessian_reset: 0
##  zero_gradients: 0
##  final_zero_gradients: 0
##  estimate_near_boundary: 0
##  s_matrix_singular: 0
##  significant_digits: 3.6
##  nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
## 
## Completed
## -- Tue Aug 16 17:21:55 2016

Load previous results NM.FOCEI <- LoadSOObject(“NONMEM_FOCEI/UseCase1_1_FOCEI.SO.xml”)

myTaskObj.FOCEI
## An object of class "taskObj"
## Slot "ESTIMATE":
## $algo
## [1] "focei"
## 
## $blocks
## $blocks[[1]]
## $blocks[[1]]$TARGET_SETTINGS
## $blocks[[1]]$TARGET_SETTINGS$INTER
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$MAXEVAL
## [1] "9999"
## 
## $blocks[[1]]$TARGET_SETTINGS$NSIG
## [1] "3"
## 
## $blocks[[1]]$TARGET_SETTINGS$SIGL
## [1] "10"
## 
## $blocks[[1]]$TARGET_SETTINGS$PRINT
## [1] "5"
## 
## $blocks[[1]]$TARGET_SETTINGS$NOABORT
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$NOPRIOR
## [1] "1"
## 
## $blocks[[1]]$TARGET_SETTINGS$FILE
## [1] "\"example1.ext\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_MATRIX
## [1] "\"R\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_PRINT
## [1] "\"E\""
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_UNCONDITIONAL
## [1] "true"
## 
## $blocks[[1]]$TARGET_SETTINGS$COV_SIGL
## [1] "12"
## 
## $blocks[[1]]$TARGET_SETTINGS$MSFO
## [1] "\"MSFO.msf\""
## 
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs
## $blocks[[1]]$TARGET_SETTINGS$blkAttrs$target
## [1] "\"NONMEM\""
## 
## 
## 
## $blocks[[1]]$.subtype
## [1] "BlockStmt"
## 
## 
## 
## 
## Slot "SIMULATE":
## NULL
## 
## Slot "EVALUATE":
## NULL
## 
## Slot "OPTIMISE":
## NULL
## 
## Slot "name":
## [1] "FOCEI_task"
scan(file.path(getwd(),"NONMEM_FOCEI",paste(uc,"_FOCEI.ctl",sep="")), what="character", sep="\n", strip.white=T)
##  [1] "; Script generated by the pharmML2Nmtran Converter v.0.4.0"                                                   
##  [2] "; Source \t: PharmML 0.8.1"                                                                                    
##  [3] "; Target \t: NMTRAN 7.3.0"                                                                                     
##  [4] "; Model \t: UseCase1_1_FOCEI"                                                                                  
##  [5] "; Model Name: UseCase1 variant 1"                                                                             
##  [6] "; Dated \t: Tue Aug 16 17:16:57 CEST 2016"                                                                     
##  [7] "$PROBLEM  Warfarin Pop PK base model"                                                                         
##  [8] "$INPUT  ID TIME WT=DROP AMT DVID DV MDV LOGTWT"                                                               
##  [9] "$DATA \"warfarin_conc.csv\" IGNORE=@"                                                                         
## [10] "$SUBS ADVAN13 TOL=9"                                                                                          
## [11] "$MODEL"                                                                                                       
## [12] "COMP (COMP1 DEFDOSE) \t;GUT"                                                                                   
## [13] "COMP (COMP2) \t;CENTRAL"                                                                                       
## [14] "$PK"                                                                                                          
## [15] "POP_CL = THETA(1)"                                                                                            
## [16] "POP_V = THETA(2)"                                                                                             
## [17] "POP_KA = THETA(3)"                                                                                            
## [18] "POP_TLAG = THETA(4)"                                                                                          
## [19] "RUV_PROP = THETA(5)"                                                                                          
## [20] "RUV_ADD = THETA(6)"                                                                                           
## [21] "BETA_CL_WT = THETA(7)"                                                                                        
## [22] "BETA_V_WT = THETA(8)"                                                                                         
## [23] "ETA_CL = ETA(1)"                                                                                              
## [24] "ETA_V = ETA(2)"                                                                                               
## [25] "ETA_KA = ETA(3)"                                                                                              
## [26] "ETA_TLAG = ETA(4)"                                                                                            
## [27] "PPV_TLAG = 0.0"                                                                                               
## [28] "MU_1 = LOG(POP_CL) + BETA_CL_WT * LOGTWT"                                                                     
## [29] "CL =  EXP(MU_1 + ETA(1)) ;"                                                                                   
## [30] "MU_2 = LOG(POP_V) + BETA_V_WT * LOGTWT"                                                                       
## [31] "V =  EXP(MU_2 + ETA(2)) ;"                                                                                    
## [32] "MU_3 = LOG(POP_KA)"                                                                                           
## [33] "KA =  EXP(MU_3 + ETA(3)) ;"                                                                                   
## [34] "MU_4 = LOG(POP_TLAG)"                                                                                         
## [35] "TLAG =  EXP(MU_4 + ETA(4)) ;"                                                                                 
## [36] "A_0(1) = 0"                                                                                                   
## [37] "A_0(2) = 0"                                                                                                   
## [38] "$DES"                                                                                                         
## [39] "GUT_DES = A(1)"                                                                                               
## [40] "CENTRAL_DES = A(2)"                                                                                           
## [41] "RATEIN_DES = 0"                                                                                               
## [42] "IF (T.GE.TLAG) THEN"                                                                                          
## [43] "RATEIN_DES = (GUT_DES*KA)"                                                                                    
## [44] "ENDIF"                                                                                                        
## [45] "CC_DES = (CENTRAL_DES/V)"                                                                                     
## [46] "DADT(1) = -(RATEIN_DES)"                                                                                      
## [47] "DADT(2) = (RATEIN_DES-((CL*CENTRAL_DES)/V))"                                                                  
## [48] "$ERROR"                                                                                                       
## [49] "EPS_Y = EPS(1)"                                                                                               
## [50] "GUT = A(1)"                                                                                                   
## [51] "CENTRAL = A(2)"                                                                                               
## [52] "RATEIN = 0"                                                                                                   
## [53] "IF (TIME.GE.TLAG) THEN"                                                                                       
## [54] "RATEIN = (GUT*KA)"                                                                                            
## [55] "ENDIF"                                                                                                        
## [56] "CC = (CENTRAL/V)"                                                                                             
## [57] "IPRED = CC"                                                                                                   
## [58] "W = RUV_ADD+RUV_PROP*IPRED"                                                                                   
## [59] "Y = IPRED+W*EPS_Y"                                                                                            
## [60] "IRES = DV - IPRED"                                                                                            
## [61] "IWRES = IRES/W"                                                                                               
## [62] "$THETA"                                                                                                       
## [63] "( 0.001 , 0.1 )\t;POP_CL"                                                                                      
## [64] "( 0.001 , 8.0 )\t;POP_V"                                                                                       
## [65] "( 0.001 , 0.362 )\t;POP_KA"                                                                                    
## [66] "( 0.001 , 1.0 )\t;POP_TLAG"                                                                                    
## [67] "( 0.0 , 0.1 )\t;RUV_PROP"                                                                                      
## [68] "( 1.0E-4 , 0.1 )\t;RUV_ADD"                                                                                    
## [69] "(0.75  FIX )\t;BETA_CL_WT"                                                                                     
## [70] "(1.0  FIX )\t;BETA_V_WT"                                                                                       
## [71] "$OMEGA BLOCK(2) CORRELATION SD"                                                                               
## [72] "(0.1 )\t;PPV_CL"                                                                                               
## [73] "(0.01 )\t;CORR_CL_V"                                                                                           
## [74] "(0.1 )\t;PPV_V"                                                                                                
## [75] "$OMEGA"                                                                                                       
## [76] "(0.1 SD )\t;PPV_KA"                                                                                            
## [77] "(0.1  FIX SD )\t;PPV_TLAG"                                                                                     
## [78] "$SIGMA"                                                                                                       
## [79] "1.0 FIX"                                                                                                      
## [80] "$EST METHOD=COND  MSFO=MSFO.msf SIGL=10 NOPRIOR=1 INTER NOABORT FILE=example1.ext MAXEVAL=9999 NSIG=3 PRINT=5"
## [81] "$COV  PRINT=E UNCONDITIONAL SIGL=12 MATRIX=R"                                                                 
## [82] "$TABLE  ID TIME AMT DVID MDV LOGTWT PRED RES WRES DV IPRED IRES IWRES Y NOAPPEND NOPRINT FILE=sdtab"          
## [83] "$TABLE  ID CL V KA TLAG ETA_CL ETA_V ETA_KA ETA_TLAG NOAPPEND NOPRINT FILE=patab"                             
## [84] "$TABLE  ID LOGTWT NOAPPEND NOPRINT FILE=cotab"

Covariance step

getPopulationParameters(NM, what="precisions") 
## $MLE
##     Parameter      MLE        SE        RSE
## 1  BETA_CL_WT 0.750000 0.0000000   0.000000
## 2   BETA_V_WT 1.000000 0.0000000   0.000000
## 3   CORR_CL_V 0.231315 0.4028760 174.167693
## 4      POP_CL 0.133780 0.0131786   9.850949
## 5      POP_KA 1.747960 1.0126300  57.932104
## 6    POP_TLAG 0.946880 0.0331714   3.503232
## 7       POP_V 8.073680 0.3690620   4.571174
## 8      PPV_CL 0.267678 0.0786570  29.384933
## 9      PPV_KA 1.016780 0.6037250  59.376168
## 10   PPV_TLAG 0.100000 0.0000000   0.000000
## 11      PPV_V 0.139330 0.0330129  23.694036
## 12    RUV_ADD 0.192779 0.0844603  43.811982
## 13   RUV_PROP 0.070990 0.0111848  15.755459