campsisnca

R-CMD-check Codecov test coverage

Analyse your simulation output using non-compartmental analysis.

Installation

Install the latest stable release as follows:

devtools::install_github("Calvagone/campsisnca")

Basic use

First import the campsisnca and gtsummary packages as follows:

library(campsisnca)
library(gtsummary)
library(gt)

Example 1: PK metrics at Day 1 and Day 7

Assume some results were simulated with Campsis (see campsis dataframe) :

campsis <- campsisnca::pk_bolus_md
campsis
## # A tibble: 5,000 × 16
##       ID  TIME   ARM   A_DEPOT A_CENTRAL A_PERIPHERAL A_OUTPUT    BW    CL    V2
##    <dbl> <dbl> <dbl>     <dbl>     <dbl>        <dbl>    <dbl> <dbl> <dbl> <dbl>
##  1     1     0     0   0              0           0        0    93.6  5.89  94.8
##  2     1     1     0   4.31e+2      538.         11.3     19.4  93.6  5.89  94.8
##  3     1     2     0   1.86e+2      722.         32.4     59.7  93.6  5.89  94.8
##  4     1     4     0   3.44e+1      741.         71.7    153.   93.6  5.89  94.8
##  5     1     6     0   6.39e+0      657.         96.3    240.   93.6  5.89  94.8
##  6     1     8     0   1.19e+0      574.        108.     316.   93.6  5.89  94.8
##  7     1    12     0   4.09e-2      447.        110.     442.   93.6  5.89  94.8
##  8     1    16     0   1.41e-3      358.         99.8    542.   93.6  5.89  94.8
##  9     1    24     0   1.67e-6      239.         72.7    688.   93.6  5.89  94.8
## 10     1    48     0   1.67e-6      316.         96.4   1588.   93.6  5.89  94.8
## # ℹ 4,990 more rows
## # ℹ 6 more variables: Q <dbl>, V3 <dbl>, KA <dbl>, CP <dbl>, OBS_CP <dbl>,
## #   Y <dbl>

Let’s define our PK metrics at Day 1 and Day 7 as follows:

# Day 1
nca_d1 <- NCAAnalysis(name = "Day 1", window = TimeWindow(0, 24), variable = "Y") %>%
  add(c(AUC(unit = "ng/mL*h"), Cmax(unit = "ng/mL"), Tmax(unit = "h"), Ctrough(unit = "ng/mL")))

# Day 7
nca_d7 <- NCAAnalysis(name = "Day 7", window = TimeWindow(144, 168), variable = "Y") %>%
  add(c(AUC(), Cmax(), Tmax(), Ctrough()))

These 2 metrics may be imported into a metrics table object, as follows. Use the calculate method to calculate the metrics in the table.

table <- NCATable() %>%
  add(c(nca_d1, nca_d7)) %>%
  calculate(campsis)

This table can be exported:

  1. To a dataframe using the export function:
table %>% export(dest = "dataframe")
## # A tibble: 24 × 4
##    metric  stat    value analysis
##    <chr>   <chr>   <dbl> <chr>   
##  1 AUC     median 134.   Day 1   
##  2 AUC     p5     102.   Day 1   
##  3 AUC     p95    168.   Day 1   
##  4 Cmax    median  10.2  Day 1   
##  5 Cmax    p5       7.85 Day 1   
##  6 Cmax    p95     13.1  Day 1   
##  7 tmax    median   2    Day 1   
##  8 tmax    p5       1    Day 1   
##  9 tmax    p95      6    Day 1   
## 10 Ctrough median   2.89 Day 1   
## # ℹ 14 more rows

When type is not specified, default value is summary. Argument type can also be summary_wide or summary_pretty. In the latter case, summary statistics are exported according to the arguments stat_display and digits provided for each metric.

  1. To an HTML table using gt:
table %>% export(dest = "gt") %>% as_raw_html()
Metric Day 1
N = 2001 		Day 7
N = 2001
AUC (ng/mL*h) 134 (102–168) 199 (131–297)
Cmax (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Ctrough (ng/mL) 2.89 (1.43–4.50) 4.18 (1.99–8.75)
1 Median (5% Centile–95% Centile)

Please note the individual metrics can also be exported to a dataframe using the export function as follows:

table %>% export(dest = "dataframe", type = "individual_wide")
## # A tibble: 400 × 6
##       id analysis   AUC  Cmax  tmax Ctrough
##    <dbl> <chr>    <dbl> <dbl> <dbl>   <dbl>
##  1     1 Day 1    122.   9.78     4    2.75
##  2     2 Day 1     92.6  6.11     6    1.70
##  3     3 Day 1    119.  10.9      2    1.97
##  4     4 Day 1    131.   9.72     2    1.99
##  5     5 Day 1    166.  10.1      4    3.79
##  6     6 Day 1    152.  10.1      2    3.04
##  7     7 Day 1    134.   9.28     2    2.70
##  8     8 Day 1    163.  13.7      2    2.91
##  9     9 Day 1    111.   9.62     2    1.94
## 10    10 Day 1    118.   9.91     4    2.90
## # ℹ 390 more rows

Example 2: PK metrics at Day 1 and Day 7 for different body weight ranges

library(dplyr)
campsis_ <- campsis %>%
  mutate(Scenario = ifelse(BW >= 75, ">=75kg patients", "<75kg patients"))

day1 <- NCAAnalysis(name = "Day 1", window = TimeWindow(0, 24), variable = "Y", strata = c(Scenario = "all")) %>%
  add(c(AUC(unit = "ng/mL*h"), Cmax(unit = "ng/mL"), Tmax(unit = "h"), Ctrough(unit = "ng/mL")))

day7 <- NCAAnalysis(name = "Day 7", window = TimeWindow(144, 168), variable = "Y", strata = c(Scenario = "all")) %>%
  add(c(AUC(), Cmax(), Tmax(), Ctrough()))

table <- NCATable() %>%
  add(c(day1, day7)) %>%
  calculate(campsis_)

table %>% export(dest = "gt") %>% as_raw_html()
Metric >=75kg patients
N = 971 		<75kg patients
N = 1031
Day 1
AUC (ng/mL*h) 122 (98–154) 144 (120–174)
Cmax (ng/mL) 10.17 (7.47–12.80) 10.26 (8.23–13.42)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Ctrough (ng/mL) 2.24 (1.25–3.90) 3.46 (2.20–4.60)
Day 7
AUC (ng/mL*h) 168 (121–236) 233 (164–305)
Cmax (ng/mL) 13.4 (9.9–18.4) 16.4 (12.6–21.1)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Ctrough (ng/mL) 3.23 (1.27–5.85) 5.95 (3.03–9.21)
1 Median (5% Centile–95% Centile) 
# Alternatively, first stratification variable can be seen in columns (use of 'tbl_merge' within gtsummary)
table <- NCATable(combine_with = "tbl_merge") %>%
  add(c(day1, day7)) %>%
  calculate(campsis_)

table %>% export(dest = "gt") %>% as_raw_html()
Metric
Day 1
Day 7
>=75kg patients
N = 971 		<75kg patients
N = 1031 		>=75kg patients
N = 971 		<75kg patients
N = 1031
AUC (ng/mL*h) 122 (98–154) 144 (120–174) 168 (121–236) 233 (164–305)
Cmax (ng/mL) 10.17 (7.47–12.80) 10.26 (8.23–13.42) 13.4 (9.9–18.4) 16.4 (12.6–21.1)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Ctrough (ng/mL) 2.24 (1.25–3.90) 3.46 (2.20–4.60) 3.23 (1.27–5.85) 5.95 (3.03–9.21)
1 Median (5% Centile–95% Centile)

Example 3: Calculate 2-compartment half-life metrics

nca <- NCAAnalysis(variable = "Y") %>%
  add(c(Thalf.2cpt.dist(), Thalf.2cpt.eff(), Thalf.2cpt.z()))

table <- NCATable() %>%
  add(nca) %>%
  calculate(campsis %>% mutate(DOSE = 1000, TAU = 24))

table %>% export(dest = "gt") %>% as_raw_html()
Metric N = 2001
t½,dist 2.57 (1.91–3.51)
t½,eff 13.4 (8.3–21.1)
t½,z 14.6 (9.6–22.7)
1 Median (5% Centile–95% Centile)

Example 4: Compute terminal half-live based on data

nca <- NCAAnalysis(variable = "Y") %>%
  add(c(Thalf(window = TimeWindow(7 * 24, 10 * 24))))

table <- NCATable() %>%
  add(nca) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric N = 2001
t½ 14.5 (9.7–22.1)
1 Median (5% Centile–95% Centile)

Example 5: Round your PK metrics

# Day 1
nca_d1 <- NCAAnalysis(name = "Day 1", window = TimeWindow(0, 24), variable = "Y") %>%
  add(AUC(digits = ~ style_sigfig(.x, 2), name = "AUC1")) %>% # At least 2 significant figures (default in gtsummary)
  add(AUC(digits = c(1, 2, 2), name = "AUC2")) %>% # Respectively 1/2/2 digit(s) after decimal for med, p5 and p95
  add(AUC(digits = ~ signif(.x, 2), name = "AUC3")) %>% # 2 significant digits only
  add(AUC(digits = list(~ round(.x / 5) * 5, ~ round(.x, 1), ~ style_number(.x)), name = "AUC4")) # 1 specific function for med, p5 and p95

# Day 7
nca_d7 <- NCAAnalysis(name = "Day 7", window = TimeWindow(144, 168), variable = "Y") %>%
  add(AUC(name = "AUC1")) %>%
  add(AUC(name = "AUC2")) %>%
  add(AUC(name = "AUC3")) %>%
  add(AUC(name = "AUC4"))

table <- NCATable()
table <- table %>%
  add(c(nca_d1, nca_d7)) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric Day 1
N = 2001 		Day 7
N = 2001
AUC1 134 (102–168) 199 (131–297)
AUC2 134.0 (102.03–167.51) 199.4 (130.64–296.69)
AUC3 130 (100–170) 200 (130–300)
AUC4 135 (102–168) 200 (130.6–297)
1 Median (5% Centile–95% Centile)

Example 6: Export custom metrics (including categorical data)

# Compute Cmax yourself using campsisnca
custom1 <- CustomMetric(fun = ~ Cmax() %>% i_value(.x, .y), name = "C_{max} custom", unit = "ng/mL")

# Check if Cmax if higher than 12 ng/mL
custom2 <- CustomMetric(
  fun = ~ (Cmax() %>% i_value(.x, .y)) > 12,
  name = "C_{max} > 12",
  unit = "%",
  categorical = TRUE
)

# Shortcut notation is also accepted
custom3 <- CustomMetric(fun = ~ Cmax > 13, name = "C_{max}> 13", unit = "%", categorical = TRUE)


# Day 1
nca_d1 <- NCAAnalysis(name = "Day 1", window = TimeWindow(0, 24), variable = "Y") %>%
  add(c(Cmax(unit = "ng/mL"), Tmax(unit = "h"), custom1, custom2, custom3))

# Day 7
nca_d7 <- NCAAnalysis(name = "Day 7", window = TimeWindow(144, 168), variable = "Y") %>%
  add(c(Cmax(), Tmax(), custom1, custom2, custom3))

table <- NCATable()
table <- table %>%
  add(c(nca_d1, nca_d7)) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric Day 1
N = 2001 		Day 7
N = 2001
Cmax (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Cmax custom (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
Cmax > 12 (%) 31 / 200 (16%) 164 / 200 (82%)
Cmax> 13 (%) 12 / 200 (6.0%) 149 / 200 (75%)
1 Median (5% Centile–95% Centile); n / N (%) 
# Alternatively, all dichotomous levels can be shown as well:
table <- NCATable(show_all_levels = TRUE)
table <- table %>%
  add(c(nca_d1, nca_d7)) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric Day 1
N = 2001 		Day 7
N = 2001
Cmax (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
tmax (h) 2.00 (1.00–6.00) 2.00 (1.00–6.00)
Cmax custom (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
Cmax > 12 (%)

    FALSE 169 / 200 (85%) 36 / 200 (18%)
    TRUE 31 / 200 (16%) 164 / 200 (82%)
Cmax> 13 (%)

    FALSE 188 / 200 (94%) 51 / 200 (26%)
    TRUE 12 / 200 (6.0%) 149 / 200 (75%)
1 Median (5% Centile–95% Centile); n / N (%)

Example 7: Geometric mean / Geometric CV

nca <- NCAAnalysis(variable = "Y") %>%
  add(c(
    AUC(unit = "ng/mL*h", stat_display = "{geomean} ({geocv}%)"),
    Cavg(unit = "ng/mL", stat_display = "{geomean} ({geocv}%)")
  ))

table <- NCATable() %>%
  add(nca) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric N = 2001
AUC (ng/mL*h) 909 (35%)
Cavg (ng/mL) 3.79 (35.17%)
1 geomean (geocv%)

Example 8: Statistics on categorical data with more than 2 levels

get_category <- function(.x, .y) {
  values <- Cmax() %>% i_value(.x, .y)
  retValue <- dplyr::case_when(
    values < 10 ~ "(1) < 10 ng/mL",
    values >= 10 & values <= 15 ~ "(2) 10-15 ng/mL",
    values > 15 ~ "(3) > 15 ng/mL",
  )
  return(retValue)
}

# Or equivalently, the 1-line purrr-style lambda expression
# get_category <- ~case_when(Cmax < 10 ~ "(1) < 10 ng/mL", Cmax >= 10 & Cmax <= 15 ~ "(2) 10-15 ng/mL", Cmax > 15 ~ "(3) > 15 ng/mL")

# Day 1
nca_d1 <- NCAAnalysis(name = "Day 1", window = TimeWindow(0, 1, time_unit = "day"), variable = "Y") %>%
  add(Cmax(unit = "ng/mL")) %>%
  add(CustomMetric(fun = get_category, name = "C_{max} categories", unit = "%", categorical = TRUE))

# Day 7
nca_d7 <- NCAAnalysis(name = "Day 7", window = TimeWindow(6, 7, time_unit = "day"), variable = "Y") %>%
  add(Cmax()) %>%
  add(CustomMetric(fun = get_category, name = "C_{max} categories", unit = "%", categorical = TRUE))

table <- NCATable()
table <- table %>%
  add(c(nca_d1, nca_d7)) %>%
  calculate(campsis)

table %>% export(dest = "gt") %>% as_raw_html()
Metric Day 1
N = 2001 		Day 7
N = 2001
Cmax (ng/mL) 10.2 (7.8–13.1) 14.8 (10.4–20.6)
Cmax categories (%)

    (1) < 10 ng/mL 87 / 200 (44%) 6 / 200 (3.0%)
    (2) 10-15 ng/mL 111 / 200 (56%) 96 / 200 (48%)
    (3) > 15 ng/mL 2 / 200 (1.0%) 98 / 200 (49%)
1 Median (5% Centile–95% Centile); n / N (%)

Example 9: Time above or below a certain threshold

In the example below, we look at the individual time above (or below) 10 ng/mL at Day 1 for the 10 first subjects.

day1 <- campsis %>%
  timerange(0, 24) %>%
  filter(ID %in% (1:10))

campsis::spaghettiPlot(day1, "Y") +
  ggplot2::geom_hline(yintercept = 10, linetype = "dashed", color = "red")

nca <- NCAAnalysis(window = TimeWindow(0, 24), variable = "Y") %>%
  add(Cmax(unit = "ng/mL*h", stat_display = "{mean}")) %>%
  add(TimeAboveLimit(limit = 10, unit = "h", stat_display = "{mean}")) %>%
  add(TimeBelowLimit(limit = 10, unit = "h", stat_display = "{mean}"))

table <- NCATable() %>%
  add(nca) %>%
  calculate(campsis %>% filter(ID %in% (1:10)))

table %>% export(dest = "dataframe", type = "individual_wide")
## # A tibble: 10 × 4
##       id  Cmax `Time above 10` `Time below 10`
##    <dbl> <dbl>           <dbl>           <dbl>
##  1     1  9.78           0                24  
##  2     2  6.11           0                24  
##  3     3 10.9            0.829            23.2
##  4     4  9.72           0                24  
##  5     5 10.1            0.427            23.6
##  6     6 10.1            0.186            23.8
##  7     7  9.28           0                24  
##  8     8 13.7            4.94             19.1
##  9     9  9.62           0                24  
## 10    10  9.91           0                24

Summary statistics can also be exported:

table %>% export(dest = "gt") %>% as_raw_html()
Metric N = 101
Cmax (ng/mL*h) 9.93
Time above 10 (h) 0.64
Time below 10 (h) 23.36
1 Mean

Example 10: Summary statistics across arms and scenarios

Let’s generate some fictitious results where two treatment arms are simulated (1g QD and 2g QD) and where scenarios are included too (Base scenario and Lower clearance).

library(campsis)
arm1 <- Arm(subjects = 24, label = "1g QD") %>%
  add(Bolus(time = 0, amount = 1000, compartment = "ABS", ii = 24, addl = 6)) %>%
  add(Observations(seq(0, 14 * 24, by = 0.1))) # 2-weeks observations

arm2 <- Arm(subjects = 24, label = "0.5g BID") %>%
  add(Bolus(time = 0, amount = 500, compartment = "ABS", ii = 12, addl = 13)) %>%
  add(Observations(seq(0, 14 * 24, by = 0.1))) # 2-weeks observations

dataset <- Dataset() %>%
  add(c(arm1, arm2))

scenario1 <- Scenario(name = "Base scenario", model = ~.x, dataset = ~.x)
scenario2 <- Scenario(
  name = "Lower clearance",
  model = ~ .x %>%
    replace(Theta(name = "CL", value = 2)),
  dataset = ~.x
)
scenarios <- Scenarios() %>% add(c(scenario1, scenario2))
results <- simulate(
  model = model_suite$pk$`2cpt_fo`,
  dataset = dataset,
  seed = 1,
  dest = "mrgsolve",
  scenarios = scenarios
)

shadedPlot(results, "CONC", colour = "ARM", strat_extra = "SCENARIO") +
  ggplot2::facet_wrap(~SCENARIO) +
  ggplot2::xlab("Time (h)") +
  ggplot2::ylab("Concentrations (ng/mL)") +
  ggplot2::labs(colour = "Arm", fill = "Arm")

NCA summary statistics are automatically calculated across all strata levels in ARM and SCENARIO.

nca <- NCAAnalysis(name = "Day 7", window = TimeWindow(144, 168), variable = "CONC") %>%
  add(AUC(unit = "ng/mL*h")) %>%
  add(Cmax(unit = "ng/mL")) %>%
  add(CustomMetric(
    fun = ~ (Cmax() %>% i_value(.x, .y)) > 30,
    name = "C_{max} > 30",
    unit = "%",
    categorical = TRUE
  )) %>%
  add(Tmax(unit = "h", digits = 2)) %>%
  add(Ctrough(unit = "ng/mL")) %>%
  add(Thalf(unit = "h", window = TimeWindow(200, "last"))) # Thalf will be estimated by dosing a linear regression on the range [200, 'last']

table <- NCATable() %>%
  add(nca) %>%
  calculate(results)

table %>% export(dest = "gt") %>% as_raw_html()
Metric 1g QD
N = 241 		0.5g BID
N = 241
Base scenario
AUC (ng/mL*h) 305 (230–543) 338 (254–467)
Cmax (ng/mL) 31 (23–53) 21 (16–25)
Cmax > 30 (%) 15 / 24 (63%) 0 / 24 (0%)
tmax (h) 0.80 (0.50–1.70) 12.80 (12.60–13.20)
Ctrough (ng/mL) 5.1 (3.1–13.2) 9.7 (6.8–15.2)
t½ (h) 14.24 (12.89–19.87) 13.44 (12.10–18.15)
Lower clearance
AUC (ng/mL*h) 457 (344–805) 506 (381–691)
Cmax (ng/mL) 37 (28–64) 28 (22–34)
Cmax > 30 (%) 21 / 24 (88%) 6 / 24 (25%)
tmax (h) 0.80 (0.50–1.90) 12.85 (12.70–13.30)
Ctrough (ng/mL) 10.4 (7.0–23.5) 16.5 (11.9–24.6)
t½ (h) 15.9 (12.5–28.5) 17.9 (12.0–26.7)
1 Median (5% Centile–95% Centile); n / N (%)

In the previous example, statistics are summarized on Day 7. In you wish specifics statistics for each one of your arms, you could also proceed sightly differently by creating 2 analyses (1 for each strata) and refer to specific arms by overriding the default strata.

nca_arm1 <- NCAAnalysis(
  name = "Last dose in '1g QD' arm",
  window = TimeWindow(144, 168),
  variable = "CONC",
  strata = c(ARM = "1g QD", SCENARIO = "all")
) %>%
  add(AUC(unit = "ng/mL*h")) %>%
  add(Cmax(unit = "ng/mL")) %>%
  add(CustomMetric(
    fun = ~ (Cmax() %>% i_value(.x, .y)) > 30,
    name = "C_{max} > 30",
    unit = "%",
    categorical = TRUE
  )) %>%
  add(Tmax(unit = "h", digits = 2)) %>%
  add(Ctrough(unit = "ng/mL")) %>%
  add(Thalf(unit = "h", window = TimeWindow(200, "last")))

nca_arm2 <- NCAAnalysis(
  name = "Last dose in '0.5 BID' arm",
  window = TimeWindow(156, 168),
  variable = "CONC",
  strata = c(ARM = "0.5g BID", SCENARIO = "all")
) %>%
  add(AUC(unit = "ng/mL*h")) %>%
  add(Cmax(unit = "ng/mL")) %>%
  add(CustomMetric(
    fun = ~ (Cmax() %>% i_value(.x, .y)) > 30,
    name = "C_{max} > 30",
    unit = "%",
    categorical = TRUE
  )) %>%
  add(Tmax(unit = "h", digits = 2)) %>%
  add(Ctrough(unit = "ng/mL")) %>%
  add(Thalf(unit = "h", window = TimeWindow(200, "last")))

table <- NCATable() %>%
  add(nca_arm1) %>%
  add(nca_arm2) %>%
  calculate(results)

table %>% export(dest = "gt") %>% as_raw_html()
Metric Base scenario
N = 241 		Lower clearance
N = 241
Last dose in ‘1g QD’ arm
AUC (ng/mL*h) 305 (230–543) 457 (344–805)
Cmax (ng/mL) 31 (23–53) 37 (28–64)
Cmax > 30 (%) 15 / 24 (63%) 21 / 24 (88%)
tmax (h) 0.80 (0.50–1.70) 0.80 (0.50–1.90)
Ctrough (ng/mL) 5.1 (3.1–13.2) 10.4 (7.0–23.5)
t½ (h) 14.2 (12.9–19.9) 15.9 (12.5–28.5)
Last dose in ‘0.5 BID’ arm
AUC (ng/mL*h) 169 (127–234) 253 (191–346)
Cmax (ng/mL) 20.7 (15.8–24.7) 27.7 (21.5–34.0)
Cmax > 30 (%) 0 / 24 (0%) 6 / 24 (25%)
tmax (h) 0.80 (0.60–1.20) 0.85 (0.70–1.30)
Ctrough (ng/mL) 9.7 (6.8–15.2) 16.5 (11.9–24.6)
t½ (h) 13.4 (12.1–18.2) 17.9 (12.0–26.7)
1 Median (5% Centile–95% Centile); n / N (%)