covariate model development
In NMproject: Script Based 'NONMEM' Model Development

This document was created by r Sys.info()['user'] at r Sys.time() in r getwd() using r R.version.string.

## DO NOT MODIFY THIS BLOCK (unless you know what you're doing)
library(knitr)
library(rprojroot)
opts_knit$set(root.dir = find_root(has_file(".Rprofile") | 
                                     is_rstudio_project | 
                                     is_r_package | 
                                     is_git_root))
opts_chunk$set(echo = TRUE)
opts_chunk$set(message = TRUE)

## LOAD PACKAGES HERE

library(NMproject)
library(dplyr)
library(ggplot2)

devtools::load_all()

Base model fit

m1 <- readRDS("Results/m1.RDS")

c1 <- m1 %>% child(run_id = "c1") %>%
  data_path("DerivedData/THEOPcov1.csv") %>%
  ## refill $INPUT because THEOPcov1.csv has different columns
  fill_input() %>%
  run_nm()

c1 %>% saveRDS("Results/base_model_pk.RDS")

Define testing relationships

# Test all parameters against a random uncorrelated covariate
dtest <- test_relations(param = c("KA", "K", "V"),
                        cov = "RND1",
                        state = "linear", 
                        continuous = TRUE) %>%
  # test LIN1 against KA
  test_relations(param = "KA",
                 cov = "LIN1", 
                 state = c("linear", "power"), 
                 continuous = TRUE) %>%
  # test LIN2 against K
  test_relations(param = "K",
                 cov = "LIN2", 
                 state = c("linear", "power"), 
                 continuous = TRUE) %>%
  # test LIN3 against V
  test_relations(param = "V",
                 cov = "LIN3", 
                 state = c("linear", "power"), 
                 continuous = TRUE) %>%
  # test BN1 against KA
  test_relations(param = "KA",
                 cov = "BN1",
                 state = "linear",
                 continuous = FALSE)

dtest

1st step forward

dsc1 <- c1 %>% covariate_step_tibble(
  run_id = "c1_f1",
  dtest = dtest,
  direction = "forward"
)

dsc1$m <- dsc1$m %>% run_nm(threads = 1)

wait_finish(dsc1$m)

dsc1 <- dsc1 %>% bind_covariate_results()

dsc1 %>% select(-m, -location)

# # uncomment to do diagnostics of top 3
# dsc1$m[1:3] %>% nm_render("Scripts/basic_gof.Rmd")

c1_f1 <- dsc1$m[1]
stopifnot(run_id(c1_f1) == "c1_f1_KA_LIN1_linear")

Decision: V - LIN3 -linear - most significant - decent gofs

2nd step forward

dsc2 <- c1_f1 %>% covariate_step_tibble(
  run_id = "c1_f2",
  dtest = dtest,
  direction = "forward"
)

dsc2$m <- dsc2$m %>% run_nm(threads = 1)

#dsc2$m <- dsc2$m %>% run_nm()

wait_finish(dsc2$m)

dsc2 <- dsc2 %>% bind_covariate_results()

dsc2 %>% select(-m, -location)

# # uncomment to do diagnostics of top 3
# dsc2$m[1:3] %>% nm_render("Scripts/basic_gof.Rmd")

c1_f2 <- dsc2$m[1]

stopifnot(run_id(c1_f2) == "c1_f2_K_LIN2_linear")

decision: KA LIN1 power - most significant even though cov step failed

1st back ward step

dsc3 <- c1_f2 %>% covariate_step_tibble(
  run_id = "c1_f2_b1",
  dtest = dtest,
  direction = "backward"
)

dsc3$m <- dsc3$m %>% run_nm(threads = 1)

wait_finish(dsc3$m)

dsc3 <- dsc3 %>% bind_covariate_results()

dsc3 %>% select(-m, -location)

Both covariates still significant. End of covariate testing

final_model <- c1_f2

summary_long(final_model, parameters = "all")

# # uncomment to get gofs of final model
# final_model %>% nm_render("Scripts/basic_gof.Rmd")

saveRDS(final_model, "Results/final_model_pk.RDS")

cov_scenarios <- final_model %>% 
  input_data() %>% 
  select(ID, LIN1, LIN2) %>%
  unique() %>%
  tidyr::gather(key = "cov", value = "value", LIN1:LIN2) %>%
  group_by(cov) %>%
  summarise(low = quantile(value, probs = 0.05),
            mid = quantile(value, probs = 0.50),
            upp = quantile(value, probs = 0.95)) %>%
  tidyr::gather(key = "type", value = "value", low:upp) %>%
  mutate(text = paste(cov, type, sep = "_"))


dplot <- final_model %>% 
  cov_forest_data(cov_scenarios)

dplot %>% cov_forest_plot()

final_model %>% param_cov_diag(param = "KA", cov = "LIN1", categorical = FALSE)

final_model %>% param_cov_diag(param = "V", cov = "LIN2", categorical = FALSE)

final_model %>% param_cov_diag(param = "V", cov = "LIN2", categorical = FALSE)

# simulation so we can do a ppc

final_model_sim <- final_model %>% child(run_id = "m1s") %>%
  update_parameters(final_model) %>%
  convert_to_simulation(subpr = 10) %>%
  run_nm()

EXPstat <- function(d, ...){ ## example statistic function
  ## arg = nonmem dataset data.frame 
  ## return data.frame with statistic column

  d <- d %>% mutate(
    LIN1_bin = Hmisc::cut2(LIN1, g = 3),
    LIN3_bin = Hmisc::cut2(LIN3, g = 3)
  )

  levels(d$LIN1_bin) <- c("LIN1_low", "LIN1_mid", "LIN1_upp")
  levels(d$LIN3_bin) <- c("LIN3_low", "LIN3_mid", "LIN3_upp")

  d %>%
    group_by(ID, LIN1_bin, LIN3_bin) %>% filter(is.na(AMT)) %>%
    summarise(
      AUC = AUC(time = TIME, conc = DV),
      CMAX = max(DV, na.rm = TRUE),
      TMAX = TIME[which.max(DV)]
    ) %>% 
    group_by(LIN1_bin, LIN3_bin) %>%
    summarise_at(c("AUC", "CMAX", "TMAX"), median, na.rm = TRUE) %>%
    tidyr::gather(AUC:TMAX, key = "type", value = "statistic")
}

Basic whisker plot - with 2 factors

dppc <- final_model_sim %>% ppc_data(EXPstat)
dppc %>% ppc_whisker_plot(type, LIN1_bin, LIN3_bin)

Only AUCs from ppc_data output and added layers to the plot

dppc %>% filter(type %in% "AUC") %>%
  ppc_whisker_plot(LIN1_bin, LIN3_bin) +
  ylab("AUC") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

With multiple factors is sometimes necessary to hold some fixed. The following set LIN3 to be it's median prior to statistic calculation.

pre_proc <- function(d){
  d$LIN3 <- median(d$LIN3)
  d
}

dppc <- final_model_sim %>% ppc_data(
  EXPstat, pre_proc = pre_proc
)

dppc %>% filter(type %in% "AUC") %>%
  ppc_histogram_plot(LIN1_bin) +
  xlab("AUC")