Prior example332_1 test

This example script is intended to to illustrate Bayesian analysis with BUGS based on example332_1.mdl following prior specification in the BayesianModels_v1.1.doc section 3.3.2.

The following steps are implemented in this workflow:

Prepared by Mike K Smith, Pfizer Ltd, Sandwich, UK 04 August 2016

Initialisation

Clear workspace and set working directory

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

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("~/UseCasesDemo/models/Priors")
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 <- "example332_1.mdl"
model <- tools::file_path_sans_ext(mdlfile)
datafile <- "Testdata.csv"

wd <- file.path(mydir,"example332v1")
dir.create(wd)

Copy the dataset and the .mdl file available under “models” into the working directory

file.copy(file.path(mydir,datafile),wd)
## [1] TRUE
file.copy(file.path(mydir,mdlfile),wd)
## [1] TRUE

Set the working directory.

setwd(wd)

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 = wd) 
cache.path <- file.path(wd, 'cache')
dir.create(cache.path)
opts_chunk$set(cache.path=cache.path)
figure.path <- file.path(wd, 'figure')
dir.create(figure.path)
opts_chunk$set(figure.path=figure.path)

List files available in working directory

list.files()
## [1] "cache"            "example332_1.mdl" "figure"          
## [4] "Testdata.csv"

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]]

Exploratory Data Analysis

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

myData <- readDataObj(myDataObj)

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

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

library(ggplot2)

Plot the data using ggplot2 library

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

plot of chunk ExploratoryAnalysis

Model Development

ESTIMATE model parameters in WINBUGS

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

# BUGS <- 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,subfolder="WinBUGS")
## -- Fri Aug 19 15:57:05 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 c20b8b77-a4f8-472e-8b13-3a938c0e0649 progress:
## Running [ ... ]
## Importing Results
## Copying the result data back to the local machine for job ID c20b8b77-a4f8-472e-8b13-3a938c0e0649...
## From C:\Users\smith_mk\AppData\Local\Temp\4\RtmpWk83qr\DDMORE.job87b4364d10c9 to C:/SEE/MDL_IDE/workspace/UseCasesDemo/models/Priors/example332v1/WinBUGS
## Done.
## 
## 
## The following main elements were parsed successfully:
##   RawResults
##   Estimation::PopulationEstimates::Bayesian
##   Estimation::PrecisionPopulationEstimates::Bayesian
##   Estimation::IndividualEstimates::Estimates
##   Estimation::PrecisionIndividualEstimates::StandardDeviation
##   Estimation::PrecisionIndividualEstimates::EstimatesDistribution
##   Estimation::PrecisionIndividualEstimates::PercentilesCI
## 
## The following MESSAGEs were raised during the job execution:
##  winbugs_message: success
## 
## Completed
## -- Fri Aug 19 15:58:10 2016
# BUGS <- LoadSOObject(file.path("WinBUGS","InstallationTest_DEQ_BUGS.SO.xml"))

getPopulationParameters(BUGS)$Bayesian
##          Parameter        Mean    Median          SDP    Perc_2.5
## 1  invOMEGA_P[1,1] 11.83153800 11.280000  4.432809729  4.94740000
## 2  invOMEGA_P[1,2]  1.74426190  1.575500  3.173392644 -4.42500000
## 3  invOMEGA_P[2,1]  1.74426190  1.575500  3.173392644 -4.42500000
## 4  invOMEGA_P[2,2] 12.20542640 11.580000  4.653772286  4.98017500
## 5        lPOP_P[1] -2.29646960 -2.296000  0.087875787 -2.48200000
## 6        lPOP_P[2]  2.08808980  2.087000  0.091539773  1.90502500
## 7     OMEGA_P[1,1]  0.10920813  0.097410  0.049376746  0.04850175
## 8     OMEGA_P[1,2] -0.01570477 -0.013285  0.034566264 -0.09218375
## 9     OMEGA_P[2,1] -0.01570477 -0.013285  0.034566264 -0.09218375
## 10    OMEGA_P[2,2]  0.10622838  0.094625  0.048818222  0.04712200
## 11         OMEGA_T  0.06266179  0.021235  0.125278138  0.00314405
## 12           POP_K  0.10100085  0.100700  0.008831009  0.08357000
## 13           POP_T  6.45692200  6.302000  1.438988570  4.04205000
## 14           POP_V  8.10324920  8.061000  0.743991578  6.72010000
## 15           TAU_T 78.43665008 47.095000 90.326820494  2.56440000
##         Perc_5     Perc_25   Perc_50      Perc_75     Perc_95   Perc_97.5
## 1   5.63215000  8.51950000 11.280000  14.48000000  19.9290000  22.0797500
## 2  -3.23290000 -0.35555000  1.575500   3.76275000   7.1525500   8.3972750
## 3  -3.23290000 -0.35555000  1.575500   3.76275000   7.1525500   8.3972750
## 4   5.72510000  8.82300000 11.580000  14.95000000  20.6600000  22.8095000
## 5  -2.44595000 -2.35200000 -2.296000  -2.23800000  -2.1510000  -2.1230000
## 6   1.93900000  2.02600000  2.087000   2.14600000   2.2440000   2.2770000
## 7   0.05381000  0.07535250  0.097410   0.13050000   0.2027000   0.2391800
## 8  -0.07310000 -0.03191750 -0.013285   0.00306425   0.0332315   0.0477965
## 9  -0.07310000 -0.03191750 -0.013285   0.00306425   0.0332315   0.0477965
## 10  0.05220900  0.07311000  0.094625   0.12640000   0.1991900   0.2294900
## 11  0.00388105  0.00929025  0.021235   0.06190000   0.2608900   0.3900375
## 12  0.08668050  0.09522000  0.100700   0.10670000   0.1163000   0.1197000
## 13  4.34800000  5.45400000  6.302000   7.33600000   8.9859000   9.5759250
## 14  6.95400000  7.58525000  8.061000   8.55200000   9.4330000   9.7430000
## 15  3.83220000 16.15250000 47.095000 107.60000000 257.6950000 318.0950000
getPopulationParameters(BUGS, what="estimates")$Bayesian # o other options,  what= precisions, intervals, all
##                       Parameter        Mean    Median
## invOMEGA_P[1,1] invOMEGA_P[1,1] 11.83153800 11.280000
## invOMEGA_P[1,2] invOMEGA_P[1,2]  1.74426190  1.575500
## invOMEGA_P[2,1] invOMEGA_P[2,1]  1.74426190  1.575500
## invOMEGA_P[2,2] invOMEGA_P[2,2] 12.20542640 11.580000
## lPOP_P[1]             lPOP_P[1] -2.29646960 -2.296000
## lPOP_P[2]             lPOP_P[2]  2.08808980  2.087000
## OMEGA_P[1,1]       OMEGA_P[1,1]  0.10920813  0.097410
## OMEGA_P[1,2]       OMEGA_P[1,2] -0.01570477 -0.013285
## OMEGA_P[2,1]       OMEGA_P[2,1] -0.01570477 -0.013285
## OMEGA_P[2,2]       OMEGA_P[2,2]  0.10622838  0.094625
## OMEGA_T                 OMEGA_T  0.06266179  0.021235
## POP_K                     POP_K  0.10100085  0.100700
## POP_T                     POP_T  6.45692200  6.302000
## POP_V                     POP_V  8.10324920  8.061000
## TAU_T                     TAU_T 78.43665008 47.095000
# it is possible also to omit $Bayesian

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

How to retrieve parameter names from the Task Properties Object

#estPars <- getTaskPropertiesObjects(mdlfile)[[1]]@ESTIMATE$blocks[[1]]$TARGET_SETTINGS$parameters
#estPars <- gsub("\\\"","", estPars)
#estPars <- gsub(" ", "", estPars)
#estPars <- unlist(strsplit(estPars, split=","))

Since there are very few parameters and they are consistent across the Priors example models:

estPars <-c("POP_K","POP_V","POP_T","OMEGA_P[1,1]","OMEGA_P[2,2]","OMEGA_P[1,2]","OMEGA_T")

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 = 1001:6000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                  Mean       SD  Naive SE Time-series SE
## POP_K         0.10100 0.008831 0.0001249      0.0003064
## POP_V         8.10325 0.743992 0.0105216      0.0267363
## POP_T         6.45692 1.438989 0.0203504      0.0555373
## OMEGA_P[1,1]  0.10921 0.049377 0.0006983      0.0008752
## OMEGA_P[2,2]  0.10623 0.048818 0.0006904      0.0008188
## OMEGA_P[1,2] -0.01570 0.034566 0.0004888      0.0005585
## OMEGA_T       0.06266 0.125278 0.0017717      0.0125051
## 
## 2. Quantiles for each variable:
## 
##                   2.5%       25%      50%      75%   97.5%
## POP_K         0.083570  0.095220  0.10070 0.106700 0.11970
## POP_V         6.723900  7.585750  8.06100 8.552000 9.74300
## POP_T         4.043950  5.454000  6.30200 7.336000 9.57307
## OMEGA_P[1,1]  0.048568  0.075357  0.09741 0.130500 0.23842
## OMEGA_P[2,2]  0.047198  0.073150  0.09463 0.126400 0.22911
## OMEGA_P[1,2] -0.091946 -0.031913 -0.01328 0.003059 0.04766
## OMEGA_T       0.003146  0.009291  0.02123 0.061900 0.38766
## trace and density plots
plot(coda_pars, ask=F)

plot of chunk MCMC_Diagnostics plot of chunk MCMC_Diagnostics

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