---
title: "Fosdagrocorat osteocalcin K-PD (Shoji 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fosdagrocorat osteocalcin K-PD (Shoji 2017)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Shoji_2017_fosdagrocorat_oc"))$reference`
* Description: `r rxode2::rxode(readModelDb("Shoji_2017_fosdagrocorat_oc"))$description`
* Article: <https://doi.org/10.1002/psp4.12201>

# Population

The osteocalcin (OSTCALC) K-PD fit used the same 321-patient phase II cohort as
the sister P1NP analysis (NCT01393639); see
[`Shoji_2017_fosdagrocorat_p1np`](Shoji_2017_fosdagrocorat_p1np.html) for the
full demographic summary. For the OSTCALC fit Shoji 2017 fixed KDE to the
P1NP-derived estimates (KDE-only re-estimation was unstable for OSTCALC data)
and fixed `Imax` to 1 for both drugs (the unconstrained OSTCALC estimate was
close to 1).

# Source trace

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lkel` (KDE fosdagrocorat) | `fixed(log(0.597))` /week | Table 2 OSTCALC: "KDE Fosdagrocorat FIX" (from P1NP fit) |
| `dlkel_pred` | `fixed(log(0.535/0.597))` | Table 2 OSTCALC: "KDE Prednisone FIX" (from P1NP fit) |
| `lkd`  | `log(0.939)` /week | Table 2 OSTCALC, "Kd" |
| `lbl`  | `log(22.2)` ng/mL | Table 2 OSTCALC, "BL" |
| `imax` | `fixed(1)` | Table 2 OSTCALC, "Imax FIX" (both drugs) |
| `ledk50` (EDK50 fos) | `log(148)` mg/week | Table 2 OSTCALC, "EDK50 Fosdagrocorat" |
| `dledk50_pred` | `log(122/148)` | Table 2 OSTCALC, "EDK50 Prednisone" |
| `hill` | `fixed(1)` | Table 2 OSTCALC, "c FIX" |
| `lrbmax` | `log(0.0276)` /mg | Table 2 OSTCALC, "RBmax" |
| `lt50`   | `log(2.24)` weeks | Table 2 OSTCALC, "T50" |
| `slp`    | `0.0675` ng/mL/week | Table 2 OSTCALC, "SLP" |
| `etalkel`   | `omega^2 = 1.5129` (123% CV) | Table 2 OSTCALC, "IIV %CV [g_KDE]" |
| `etaledk50` | `omega^2 = 0.0650` (25.5% CV) | Table 2 OSTCALC, "IIV %CV [g_EDK50]" |
| `etalbl`    | `omega^2 = 0.1901` (43.6% CV) | Table 2 OSTCALC, "IIV %CV [g_BL]" |
| `etaslp`    | `omega^2 = 0.338^2 = 0.1142` | Table 2 OSTCALC, "IIV SD [g_SLP] = 0.338" |
| `propSd`    | `0.141` | Table 2 OSTCALC, "Residual variability %CV [e] = 14.1" |
| `d/dt(depot_kpd)`  | -kel * depot_kpd | Methods, K-PD model equations |
| `d/dt(effect)` | ks * rebound * inhibition - kd * effect | Methods, K-PD model + Results rebound equation |
| `OSTCALC = effect + slp_i * t` | Observation = response + linear placebo trend | Methods, F(ij) equation |

# Virtual cohort

The cohort and dosing scheme mirror the P1NP vignette: 200 virtual subjects
per trial arm, 8 weeks active treatment followed by a tapered 4-week period.

```{r cohort}
set.seed(20170527L)

n_per_arm <- 200L

make_arm <- function(arm_label, dose_qd_mg, drug_pred, n,
                     dose_reduced_mg, id_offset) {
  active_rate  <- 7 * dose_qd_mg
  taper_a_rate <- 3.5  * dose_reduced_mg
  taper_b_rate <- 2.33 * dose_reduced_mg

  ids <- id_offset + seq_len(n)
  dose_ev <- if (dose_qd_mg > 0) {
    bind_rows(
      tibble(id = ids, time = 0,  amt = active_rate * 8,
             rate = active_rate, evid = 1L, cmt = "depot_kpd"),
      tibble(id = ids, time = 8,  amt = taper_a_rate * 2,
             rate = taper_a_rate, evid = 1L, cmt = "depot_kpd"),
      tibble(id = ids, time = 10, amt = taper_b_rate * 2,
             rate = taper_b_rate, evid = 1L, cmt = "depot_kpd")
    )
  } else {
    tibble(id = integer(), time = numeric(), amt = numeric(),
           rate = numeric(), evid = integer(), cmt = character())
  }

  obs_ev <- expand.grid(id = ids,
                        time = c(0, 2, 4, 6, 8, 10, 12, 13)) |>
    as_tibble() |>
    mutate(amt = NA_real_, rate = NA_real_, evid = 0L, cmt = NA_character_)

  bind_rows(dose_ev, obs_ev) |>
    arrange(id, time, desc(evid)) |>
    mutate(arm = arm_label, DOSE = dose_qd_mg, DRUG_PRED = drug_pred)
}

events <- bind_rows(
  make_arm("Placebo",        0,  0, n_per_arm,  0, id_offset =   0L),
  make_arm("Fos 1 mg",       1,  0, n_per_arm,  1, id_offset = 200L),
  make_arm("Fos 5 mg",       5,  0, n_per_arm,  1, id_offset = 400L),
  make_arm("Fos 10 mg",     10,  0, n_per_arm,  1, id_offset = 600L),
  make_arm("Fos 15 mg",     15,  0, n_per_arm,  1, id_offset = 800L),
  make_arm("Pred 5 mg",      5,  1, n_per_arm,  5, id_offset = 1000L),
  make_arm("Pred 10 mg",    10,  1, n_per_arm,  5, id_offset = 1200L)
)

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Shoji_2017_fosdagrocorat_oc")

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("arm", "DOSE", "DRUG_PRED")
) |> as.data.frame()

sim_typ <- rxode2::rxSolve(
  rxode2::zeroRe(mod), events = events,
  keep = c("arm", "DOSE", "DRUG_PRED")
) |> as.data.frame()
```

# Replicate Figure 2 (lower row): VPC of OSTCALC percent change from baseline

```{r figure-2-oc}
arm_order <- c("Placebo",
               "Fos 1 mg", "Fos 5 mg", "Fos 10 mg", "Fos 15 mg",
               "Pred 5 mg", "Pred 10 mg")

vpc_summary <- sim |>
  filter(time %in% c(0, 2, 4, 6, 8, 10, 12, 13)) |>
  group_by(arm, id) |>
  mutate(oc_baseline = first(OSTCALC[time == 0])) |>
  ungroup() |>
  mutate(cfb_pct = 100 * (OSTCALC - oc_baseline) / oc_baseline) |>
  group_by(arm, time) |>
  summarise(
    p10 = quantile(cfb_pct, 0.10, na.rm = TRUE),
    p50 = quantile(cfb_pct, 0.50, na.rm = TRUE),
    p90 = quantile(cfb_pct, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(arm = factor(arm, levels = arm_order))

ggplot(vpc_summary, aes(time, p50)) +
  geom_ribbon(aes(ymin = p10, ymax = p90), alpha = 0.25, fill = "darkorange") +
  geom_line(color = "darkorange") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  geom_vline(xintercept = 8, linetype = "dotted", colour = "grey50") +
  facet_wrap(~ arm, ncol = 4) +
  labs(x = "Time (weeks)",
       y = "Osteocalcin change from baseline (%)",
       title = "Replicates Figure 2 (bottom row) of Shoji 2017",
       subtitle = "Median (line) and 10-90% predictive interval (shaded) per arm",
       caption = "Vertical dotted line at week 8 marks the end of the active treatment period.")
```

# Replicate Table 3: simulated median OSTCALC %CFB at week 8

```{r table-3-oc}
published_table3 <- tibble::tribble(
  ~arm,         ~published_median_pct,
  "Placebo",                       1.8,
  "Fos 1 mg",                      0.1,
  "Fos 5 mg",                     -6.7,
  "Fos 10 mg",                   -12.6,
  "Fos 15 mg",                   -16.8,
  "Pred 5 mg",                    -9.7,
  "Pred 10 mg",                  -16.9
)

simulated_table3 <- sim |>
  group_by(id, arm) |>
  summarise(
    oc_baseline = first(OSTCALC[time == 0]),
    oc_week8    = first(OSTCALC[time == 8]),
    cfb_pct     = 100 * (oc_week8 - oc_baseline) / oc_baseline,
    .groups = "drop"
  ) |>
  group_by(arm) |>
  summarise(simulated_median_pct = median(cfb_pct, na.rm = TRUE),
            .groups = "drop")

comparison <- published_table3 |>
  left_join(simulated_table3, by = "arm") |>
  mutate(arm = factor(arm, levels = arm_order)) |>
  arrange(arm) |>
  mutate(delta = simulated_median_pct - published_median_pct)

comparison |>
  dplyr::rename(
    "Arm"                             = arm,
    "Published median %CFB (Table 3)" = published_median_pct,
    "Simulated median %CFB"           = simulated_median_pct,
    "Difference (pp)"                 = delta
  ) |>
  knitr::kable(digits = 1,
             caption = "Osteocalcin percent change from baseline at week 8 -- published vs simulated.")
```

# Typical-value parameter-recovery checks

```{r recovery-oc}
sim_typ_summary <- sim_typ |>
  filter(time %in% c(0, 8)) |>
  group_by(arm, time) |>
  summarise(OSTCALC_typ = first(OSTCALC), .groups = "drop") |>
  pivot_wider(names_from = time, values_from = OSTCALC_typ, names_prefix = "wk") |>
  mutate(cfb_pct_typ = 100 * (wk8 - wk0) / wk0,
         arm = factor(arm, levels = arm_order)) |>
  arrange(arm)

sim_typ_summary |>
  dplyr::rename(
    "Arm"               = arm,
    "OSTCALC wk 0 (typical)" = wk0,
    "OSTCALC wk 8 (typical)" = wk8,
    "%CFB typical"      = cfb_pct_typ
  ) |>
  knitr::kable(digits = 2,
             caption = "Typical-value (zero-RE) osteocalcin at week 8 per arm.")
```

Expected typical responses:

* Placebo: OSTCALC = BL + SLP * 8 = 22.2 + 0.0675 * 8 = `r 22.2 + 0.0675 * 8` ng/mL.
* Fosdagrocorat 10 mg q.d.: ~ -13% CFB at week 8 (paper Discussion / Table 3).

# Assumptions and deviations

- **Observation variable naming.** The single-output observation is named
  `OSTCALC` (the paper's name) rather than the canonical `Cc`. Same accepted
  deviation as the P1NP sister model and the wider codebase pattern for
  paper-named biomarker outputs.
- **Drug arm switching via reparameterization.** Same convention as the
  P1NP sister model (base = fosdagrocorat; `dlkel_pred` and `dledk50_pred`
  log-ratio offsets recover the published prednisone values when
  `DRUG_PRED = 1`). For OSTCALC the typical-value KDE is fixed to the P1NP
  estimates -- both `lkel` and `dlkel_pred` are wrapped in `fixed()`.
- **Imax fixed to 1 for both drugs.** The OSTCALC fit set `Imax = 1` for both
  fosdagrocorat and prednisone because the unconstrained estimate was
  close to 1 (paper Discussion). The inhibition term `1 - imax * ir^c /
  (edk50^c + ir^c)` therefore approaches 0 (full inhibition) at very high
  IR; the typical OSTCALC %CFB at week 8 for the 15 mg fosdagrocorat arm
  (-16.8% published; ~ -17% simulated) reflects this nearly-complete
  inhibition combined with the slow OSTCALC turnover (Kd = 0.939 /week, longer
  half-life than the 8-week treatment window).
- **Independent (non-block) IIV.** Shoji 2017 used independent IIVs on
  KDE, EDK50, BL for the OSTCALC fit (paper Methods: "To reduce estimation
  instability, individual random effect parameters for KDE, EDK50, and
  BL were assumed to be independent."). The model encodes this with
  three diagonal `etal*` declarations (no block). The IIV CV%s for KDE
  in OSTCALC (123%) are larger than for KDE in P1NP (95%) because the OSTCALC fit
  fits g_KDE to OSTCALC residuals alone; this is the paper's published value
  and is preserved here.
- **Taper-period dosing approximation.** Same approximation as the P1NP
  vignette (zero-order rates of 3.5 and 2.33 doses/week of the reduced
  dose for weeks 8-10 and 10-12).
