The challenge

Use data visualization to characterize the PK behavior of the drug.

Ideas

There is no correct answer, but you could explore:

The data

Simulated pharmacokinetic (PK) data:

Dataset is simple: variables for Subject, Time (in hours) post-dose, and concentration.

Spaghetti plots

ggplot(data = PK_data, aes(x = time_hr, y = concentration_ng_ml)) +
  geom_line(aes(group = subject, colour = subject)) +
  scale_colour_discrete_sequential(palette = "Red-Blue") +
  scale_x_continuous(breaks = time, labels = time) +
  labs(title = "Individual concentration-time profiles (n=50)",
       x = "Time post-dose (hours)", y = "Concentration (ng/ml)") +
  scale_y_continuous(limits = c(0, 90), breaks = seq(0, 90, by=10)) +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Mean concentration-time (linear scale) with IQR and P10-P90 bands

mean_table <- PK_data %>%
  group_by(time_hr) %>%
  summarise(Mean = mean(concentration_ng_ml),
            Median = median(concentration_ng_ml),
            SD = sd(concentration_ng_ml),
            Q25 = quantile(concentration_ng_ml, 0.25),
            Q75 = quantile(concentration_ng_ml, 0.75),
            Q10 = quantile(concentration_ng_ml, 0.10),
            Q90 = quantile(concentration_ng_ml, 0.90))

ggplot(data = mean_table, aes(x = time_hr, y = Median)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.3, fill = "skyblue") +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.3, fill = "skyblue2") +
  geom_line(aes(group = 1), colour = "darkblue") +
  geom_point(colour = "darkblue") +
  scale_x_continuous(breaks = time, labels = time) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80)) +
  labs(title = "Concentration-time, IQR, P10-P90, linear scale (n=50)",
       x = "Time post-dose (hours)", y = "Concentration (ng/ml)") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Mean concentration-time (semi-logarithmic scale)

PK_data$log_concentration <- log(PK_data$concentration_ng_ml)

mean_semi_table <- PK_data %>%
  group_by(time_hr) %>%
  summarise(LogMean = mean(log_concentration),
            LogSD = sd(log_concentration),
            Median = median(log_concentration),
            Q25 = quantile(log_concentration, 0.25),
            Q75 = quantile(log_concentration, 0.75),
            Q10 = quantile(log_concentration, 0.10),
            Q90 = quantile(log_concentration, 0.90))

ggplot(data = mean_semi_table, aes(x = time_hr, y = LogMean)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.3, fill = "skyblue") +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.3, fill = "skyblue2") +
  geom_line(aes(group = 1), colour = "darkblue") +
  geom_point(colour = "darkblue") +
  scale_x_continuous(breaks = time, labels = time) +
  # scale_y_continuous(breaks = c(0, 1, 2, 3, 4)) +
  labs(title = "Concentration-time, IQR, P10-P90, semi-logarithmic scale (n=50)",
       x = "Time post-dose (hours)", y = "Concentration (ng/ml), log scale") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Variability by timepoint (boxplots)

Spread and skew at each sampling time. The box shows the IQR; every jittered point is one subject.

PK_data$time_hr_c <- as.factor(PK_data$time_hr)

ggplot(data = PK_data, aes(x = time_hr_c, y = concentration_ng_ml)) +
  geom_point(position = position_jitter(width = 0.15), colour = "grey") +
  geom_boxplot(fill = "skyblue", alpha = 0.3) +
  labs(title = "Variability by timepoint (n=50)",
       x = "Time post-dose (hours)", y = "Concentration (ng/ml)") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Maximum concentration per participant, Tmax, Cmax vs Tmax

cmax_table <- PK_data %>%
  group_by(subject) %>%
  summarise(Cmax = max(concentration_ng_ml))

tmax_table <- merge(PK_data, cmax_table, by="subject") %>%
  filter(concentration_ng_ml == Cmax) 

tmax_table$Tmax <- tmax_table$time_hr
tmax_table$Tmax <- as.factor(tmax_table$Tmax)
tmax_table <- tmax_table %>%
  select(subject, Tmax, Cmax)

ggplot(data = tmax_table, aes(x = Tmax, y = Cmax)) +
  geom_boxplot(fill = "skyblue", alpha = 0.3) +
  geom_point(position = position_jitter(width = 0.15), colour = "#D55E00") +
  labs(title = "Individual peak concentration by peak time (n=50)",
       x = "Tmax - Time of individual peak (hours)", y = "Cmax - Individual peak concentration (ng/ml)") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Normalised profile shape C/Cmax

Normalising each profile by its subject-level peak (C/Cmax) highlights shape differences independent of magnitude.

PK_data_2 <- merge(PK_data, tmax_table, by = "subject")

PK_data_2$CCmax <- PK_data_2$concentration_ng_ml / PK_data_2$Cmax

ggplot(data = PK_data_2, aes(x = time_hr, y = CCmax)) +
  # geom_hline(aes(yintercept = 0.5), linetype = "dashed", colour = "grey") +
  geom_line(aes(group = subject, colour = subject)) +
  scale_colour_discrete_sequential(palette = "Red-Blue") +
  scale_x_continuous(breaks = time, labels = time) +
  labs(title = "Individual normalised by subject-level peak profiles (n=50)",
       x = "Time post-dose (hours)", y = "Normalised concentration (C/Cmax)") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))

Calculate AUC and apparent terminal half-life

Using example from https://vis-sig.github.io/blog/posts/2026-07-08-wonderful-wednesday-july-2026/#example2%20code

# --- Per-subject exposure: trapezoidal AUC over the observed window ---
auc_by_subject <- PK_data |>
  dplyr::group_by(subject) |>
  dplyr::arrange(time_hr, .by_group = TRUE) |>
  dplyr::summarise(
    AUC = sum(diff(time_hr) * (head(concentration_ng_ml, -1) + tail(concentration_ng_ml, -1)) / 2),
    .groups = "drop"
  )
t_first <- min(PK_data$time_hr)
t_last <- max(PK_data$time_hr)

# --- Apparent terminal half-life from the log-linear tail (time >= 8 h) ---
est_lambda <- function(t, conc) {
  ok <- is.finite(conc) & conc > 0
  if (sum(ok) < 2) return(NA_real_)
  slope <- unname(coef(stats::lm(log(conc[ok]) ~ t[ok]))[2])
  if (!is.finite(slope) || slope >= 0) return(NA_real_)
  -slope
}

half_life_by_subject <- PK_data |>
  dplyr::filter(time_hr >= 8) |>
  dplyr::group_by(subject) |>
  dplyr::summarise(lambda_z = est_lambda(time_hr, concentration_ng_ml), .groups = "drop") |>
  dplyr::mutate(t_half = log(2) / lambda_z)

Cmax vs AUC

cmax_auc_table <- merge(cmax_table, auc_by_subject, by="subject")

ggplot(data = cmax_auc_table, aes(x = AUC, y = Cmax)) +
  geom_point() +
  labs(title = "Exposure Relationship", 
       subtitle = "Max Concentration (Cmax) versus Area Under the Curve (AUC)",
       x = "AUC (ng·hours/mL)", y = "Cmax (ng/mL)") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 8))