3 Dataset Creation

library(knitr)
library(PKPDmisc)
library(tidyverse)
library(vancomycin)
library(mrgsolve)
source("../scripts/generate_demographics.R")
#> 
#> Attaching package: 'magrittr'
#> The following object is masked from 'package:purrr':
#> 
#>     set_names
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
models <- source("../models/models.R")$value
#> model dir set to C:/Users/devin/Documents/Repos/bayesiantutorial/models
#> cache location set to C:\Users\devin\Documents\Repos\bayesiantutorial\models\.modelcache
#> Compiling vanc_stockmann ...
#> done.

3.1 Generate data for mrgsolve

start with a baseline of having 100 individuals worth of data, can scale to different amounts of individuals later

set.seed(1234567)
NIDS <- 100
demogs <- generate_demographics(NIDS) %>% mutate(ID = 1:nrow(.), 
                                                 PMA = PMA/7 # in weeks for stockmann model
                                                 )
for_dosing <- demogs %>% 
                  mutate(
                      CMT = 1, 
                      EVID = 1,
                      TIME = 0
                  ) %>% select(ID, everything())
all_dosing <- for_dosing %>% group_by(ID) %>%
    nest() %>%
    mutate(dose = map(data, ~ dosing_6bin_pna_wt(.$AGE, .$WT))
           ) %>% unnest() 

Dose for 8 days

all_dosing <- all_dosing %>% 
  mutate(addl = as.integer(ceiling(24/II*8)))
kable(head(all_dosing, n = 15))
ID WT CRE AGE PMA HT SEX CMT EVID TIME DOSE II addl
1 1.74 0.972 16.13 35.7 39.5 FEMALE 1 1 0 26.1 12 16
2 2.54 0.472 31.81 40.1 53.2 FEMALE 1 1 0 38.1 8 24
3 3.62 0.377 16.79 37.0 48.1 FEMALE 1 1 0 54.3 8 24
4 1.78 0.302 8.86 33.0 38.3 FEMALE 1 1 0 26.7 12 16
5 3.34 0.254 15.48 41.4 39.6 FEMALE 1 1 0 50.1 8 24
6 2.66 0.604 32.57 35.9 51.9 FEMALE 1 1 0 39.9 8 24
7 2.32 0.674 3.53 38.1 45.8 MALE 1 1 0 34.9 12 16
8 3.31 0.375 58.36 42.0 45.5 FEMALE 1 1 0 49.6 8 24
9 2.19 0.774 6.37 30.9 28.6 MALE 1 1 0 32.9 12 16
10 2.33 1.289 12.62 35.4 44.3 FEMALE 1 1 0 34.9 8 24
11 2.52 0.219 65.33 41.2 46.6 MALE 1 1 0 37.8 8 24
12 6.92 0.558 96.26 60.0 69.0 MALE 1 1 0 103.8 8 24
13 1.47 0.932 7.76 31.4 34.1 FEMALE 1 1 0 22.1 12 16
14 3.26 0.262 25.95 38.5 60.2 FEMALE 1 1 0 48.9 8 24
15 3.14 0.391 68.77 41.7 44.7 MALE 1 1 0 47.2 8 24
vanc_stockmann <- models$use("vanc_stockmann")

Normalize column capitalizations to make mrgsolve happy

mrg_data <- bind_cols(
  all_dosing %>% select(ID),
  all_dosing %>% select(everything(), -ID) %>% lowercase_names()
  ) %>% mutate(amt = dose)
raw_sim <- vanc_stockmann %>% data_set(mrg_data) %>% mrgsim(end = 10*24, delta = 2) %>%
  as.data.frame %>% as_data_frame %>% capitalize_names()
#> Dropping non-numeric columns: sex

3.1.1 Check distribution of trough values at SS

peak_trough <- raw_sim %>% filter(between(TIME, 150, 200)) %>%
  select(ID, TIME, CP, DV) %>%
  gather(obs, value, CP, DV) %>%
  group_by(ID, obs) %>%
  summarize(ctrough = min(value),
            cmax = max(value))

Peak and trough values are representative of real observations

peak_trough %>%
  gather(sample, value, ctrough, cmax) %>%
  ggplot(aes(x = value)) + 
  geom_density(aes(fill = obs), alpha = 0.6) +
  facet_wrap(~sample, scales = "free") + theme_bw() +
  base_theme()

  • take simulated values and prepare for nonmem estimation
from_sim <- raw_sim %>% 
  distinct(ID, TIME, DV)

from_orig <-  mrg_data %>% 
  capitalize_names() %>% 
  select(
    ID,
    TIME,
    WT,
    CRE,
    AGE,
    PMA, 
    HT,
    SEX,
    CMT,
    DOSE,
    II,
    ADDL
  )
data_for_nonmem <- left_join(from_sim, from_orig) %>% 
  group_by(ID) %>%
  fill(WT:CMT) %>%
  mutate(
    EVID = ifelse(DV == 0, 1, 0)
  ) %>% 
  ## don't need to go too far after final dose for estimation purposes
  ## also don't n
  filter(TIME < 216) %>% 
  ## don't need 0 concentration time point as not realistic
  distinct(ID, TIME, .keep_all = TRUE)
#> Joining, by = c("ID", "TIME")

3.1.2 write data for nonmem

write_nonmem(data_for_nonmem, "../modeling/mdata/rich_100.csv")
devtools::session_info()
#> Session info --------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.3.2 (2016-10-31)
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  English_United States.1252  
#>  tz       America/New_York            
#>  date     2016-12-06
#> Packages ------------------------------------------------------------------
#>  package   * version date       source        
#>  backports   1.0.4   2016-10-24 CRAN (R 3.3.2)
#>  bookdown    0.2     2016-11-12 CRAN (R 3.3.2)
#>  devtools    1.12.0  2016-06-24 CRAN (R 3.3.2)
#>  digest      0.6.10  2016-08-02 CRAN (R 3.3.2)
#>  evaluate    0.10    2016-10-11 CRAN (R 3.3.2)
#>  htmltools   0.3.5   2016-03-21 CRAN (R 3.3.2)
#>  httpuv      1.3.3   2015-08-04 CRAN (R 3.3.2)
#>  knitr       1.15    2016-11-09 CRAN (R 3.3.2)
#>  magrittr    1.5     2014-11-22 CRAN (R 3.3.2)
#>  memoise     1.0.0   2016-01-29 CRAN (R 3.3.2)
#>  mime        0.5     2016-07-07 CRAN (R 3.3.2)
#>  miniUI      0.1.1   2016-01-15 CRAN (R 3.3.2)
#>  R6          2.2.0   2016-10-05 CRAN (R 3.3.2)
#>  Rcpp        0.12.8  2016-11-17 CRAN (R 3.3.2)
#>  rmarkdown   1.2     2016-11-21 CRAN (R 3.3.2)
#>  rprojroot   1.1     2016-10-29 CRAN (R 3.3.2)
#>  shiny       0.14.2  2016-11-01 CRAN (R 3.3.2)
#>  stringi     1.1.2   2016-10-01 CRAN (R 3.3.2)
#>  stringr     1.1.0   2016-08-19 CRAN (R 3.3.2)
#>  withr       1.0.2   2016-06-20 CRAN (R 3.3.2)
#>  xtable      1.8-2   2016-02-05 CRAN (R 3.3.2)
#>  yaml        2.1.13  2014-06-12 CRAN (R 3.3.2)