---
title: "Pioglitazone (Stringer 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pioglitazone (Stringer 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Stringer_2015_pioglitazone"))$reference`
* Description: `r rxode2::rxode(readModelDb("Stringer_2015_pioglitazone"))$description`
* Article: <https://doi.org/10.1089/dia.2014.0222>

Stringer et al. (2015) built a population pharmacodynamic model of fasting
plasma glucose (FPG) and glycosylated hemoglobin (HbA1c) in Japanese type 2
diabetes mellitus (T2DM) patients over a 2.5-4 year randomized study
comparing pioglitazone against a non-thiazolidinedione (non-TZD) oral
glucose-lowering drug background. The model is a simultaneous cascading
indirect-response system for FPG and HbA1c with proportional disease
progression on FPG and a linear-in-time FPG-independent input to HbA1c. No
plasma-concentration data were collected; the drug effect is expressed as a
time-driven Emax model whose Emax and ET50 differ between the two treatment
cohorts (Stringer 2015 Methods; Table 2).

# Population

The analysis pooled n = 587 Japanese T2DM patients (293 pioglitazone,
294 control) from a single multicenter, prospective, randomized,
open-label, blinded-end-point cardiovascular-outcome study
(Stringer 2015 Table 1). Median age was 58 years (range 35-74), median
body weight 68-69 kg, median BMI 26.2-26.5 kg/m^2, and 37.8% female.
Baseline HbA1c > 6.9% (NGSP) was an inclusion criterion; median baseline
FPG was 153-157 mg/dL and median baseline HbA1c was 7.6-7.9% (NGSP).
Patients received titrated treatment aimed at HbA1c < 6.9% NGSP;
concomitant non-TZD oral glucose-lowering drugs at baseline covered
sulfonylureas (73-82%), alpha-glucosidase inhibitors (36-56%),
biguanides (43-68%), and rapid-acting insulin secretagogues (7-13%).
Median treatment duration was 3.14 years (maximum 3.9 years).

The population metadata is available programmatically via
`readModelDb("Stringer_2015_pioglitazone")()$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/pharmacodynamics/Stringer_2015_pioglitazone.R`. The table
below collects the primary structural sources.

| Element                                         | Value / form              | Source location       |
|-------------------------------------------------|---------------------------|-----------------------|
| Eq. 1: `KinGDP = BSLG * KoutG * (1 + FPGDP * t)`| ODE production term       | Stringer 2015 Eq. 1   |
| Eq. 2: `d(FPG)/dt = KinGDP - KoutG * (1 + DEF) * FPG` | FPG ODE             | Stringer 2015 Eq. 2   |
| Eq. 3: `d(HbA1c)/dt = FPGind + KinH * FPG^c - KoutH * HbA1c` | HbA1c ODE  | Stringer 2015 Eq. 3   |
| Eq. 4: `FPGind = KinZ * (1 + DPind * t)`        | HbA1c FPG-independent input | Stringer 2015 Eq. 4 |
| Eq. 5: `DEF = Emax * t / (ET50 + t)`            | Time-driven Emax          | Stringer 2015 Eq. 5   |
| `BSL_FPG (females)` = 156.0 mg/dL               | fixed effect              | Stringer 2015 Table 2 |
| `KoutG` = 0.0089 days^-1                        | fixed effect              | Stringer 2015 Table 2 |
| `BSL_HbA1c` = 7.83 %                            | fixed effect              | Stringer 2015 Table 2 |
| `KoutH` = 0.072 days^-1                         | fixed effect              | Stringer 2015 Table 2 |
| `KinZ` = 0.28 days^-1 (of HbA1c %)              | fixed effect              | Stringer 2015 Table 2 |
| `c` (power on FPG for HbA1c production) = 1.91  | fixed effect              | Stringer 2015 Table 2 |
| `Emax` (pioglitazone) = 17.3%                   | fixed effect              | Stringer 2015 Table 2 |
| `Emax` (control) = 8.4%                         | fixed effect              | Stringer 2015 Table 2 |
| `ET50` (pioglitazone) = 0 days FIX              | fixed at 0                | Stringer 2015 Table 2 |
| `ET50` (control) = 49.2 days                    | fixed effect              | Stringer 2015 Table 2 |
| `FPGDP` = 0.017 year^-1                         | fixed effect              | Stringer 2015 Table 2 |
| `DPind` = 0.03 year^-1                          | fixed effect              | Stringer 2015 Table 2 |
| Sex effect on BSL_FPG = 0.05 (+5% for males)    | fixed effect              | Stringer 2015 Table 2 / Results Covariate analysis |
| IIV: OMEGA BLOCK on BSL_FPG and BSL_HbA1c       | omega^2 0.03, cov 0.01, omega^2 0.01 | Stringer 2015 Table 2 |
| IIV: `FPGDP` additive random effect             | omega^2 = 0.004           | Stringer 2015 Table 2 / Methods |
| IIV: `Emax` shared log-normal                   | omega^2 = 0.75            | Stringer 2015 Table 2 |
| Residual error FPG (proportional)               | 14.4% CV                  | Stringer 2015 Table 2 |
| Residual error HbA1c (proportional)             | 5.8% CV                   | Stringer 2015 Table 2 |

# Validation strategy

This is a pure PD turnover model with no PK data and no exogenous dose
records (Stringer 2015 Methods: "no pharmacokinetic data were collected").
PKNCA-based NCA is not the right validation target for such a model; the
established alternatives (Steady-state hold, perturbation-recovery,
mass-balance / flux check, dimensional analysis) apply here (see
`references/endogenous-validation.md`). Below we run:

1. A steady-state hold with the drug effect and both disease-progression
   terms turned off.
2. A per-cohort typical-value simulation reproducing Stringer 2015
   Figure 2 (median FPG and HbA1c time profiles for the pioglitazone and
   control groups over 5 years).
3. A tabular comparison of derived summary values against those reported
   in the paper (maximum FPG reduction, HbA1c-at-2.5-years, and the
   longer-term projected 5-year FPG differences from Stringer 2015
   Results).

## 1. Steady-state hold

At `t = 0` with no drug effect (`DEF = 0` because the Emax model's
numerator carries `t`) and both `FPGDP * t = 0` and `DPind * t = 0`, the
right-hand side of the FPG and HbA1c ODEs should evaluate to zero. In
this model KinH is derived per subject from the HbA1c steady-state
balance `KinH = (KoutH * BSL_HbA1c - KinZ) / BSL_FPG^c`, so the initial
condition is consistent by construction -- and we can check that
condition numerically by running an override where `FPGDP`, `DPind`,
`Emax_pio`, and `Emax_ctrl` are all set to zero (so no time-dependent
perturbation acts on the state) and confirming that the state stays at
its baseline.

```{r ss-check}
mod <- readModelDb("Stringer_2015_pioglitazone")
mod_typ <- rxode2::zeroRe(mod)

# Override at solve time: zero disease progression and zero drug effect,
# so the ODE right-hand side stays at zero for all t.
ev_ss <- rxode2::et(id = 1L, seq(0, 5 * 365, by = 30)) |>
  as.data.frame() |>
  mutate(SEXF = 1L, TRT = 0L)

ss <- rxode2::rxSolve(
  mod_typ, ev_ss,
  params = c(fpgdp = 0, dpind = 0,
             lemax_pio  = log(1e-12),
             lemax_ctrl = log(1e-12))
) |> as.data.frame()

stopifnot(
  isTRUE(all.equal(range(ss$fpg),   c(156, 156),   tolerance = 1e-4)),
  isTRUE(all.equal(range(ss$hba1c), c(7.83, 7.83), tolerance = 1e-4))
)
cat("FPG   range at SS:", format(range(ss$fpg),   digits = 6), "\n")
cat("HbA1c range at SS:", format(range(ss$hba1c), digits = 6), "\n")
```

## 2. Reproduce Figure 2 -- typical-value simulation

Stringer 2015 Figure 2 shows the simulated median FPG (panel a) and
HbA1c (panel b) trajectories for the two treatment groups over 5 years.
Here we run a typical-value simulation (`zeroRe`) for four cohorts
covering both treatment and sex, then reproduce panels 2a and 2b.

```{r typical-cohorts}
make_typical_cohort <- function(id, trt, sexf, times) {
  rxode2::et(id = id, times) |>
    as.data.frame() |>
    mutate(TRT = trt, SEXF = sexf)
}

times_days <- seq(0, 5 * 365, by = 15)
cohorts <- bind_rows(
  make_typical_cohort(1L, trt = 1L, sexf = 1L, times = times_days) |>
    mutate(cohort = "Pioglitazone (F)"),
  make_typical_cohort(2L, trt = 0L, sexf = 1L, times = times_days) |>
    mutate(cohort = "Control (F)"),
  make_typical_cohort(3L, trt = 1L, sexf = 0L, times = times_days) |>
    mutate(cohort = "Pioglitazone (M)"),
  make_typical_cohort(4L, trt = 0L, sexf = 0L, times = times_days) |>
    mutate(cohort = "Control (M)")
)
stopifnot(!anyDuplicated(unique(cohorts[, c("id", "time", "evid")])))

sim_typ <- rxode2::rxSolve(mod_typ, events = cohorts,
                           keep = c("cohort", "TRT", "SEXF")) |>
  as.data.frame() |>
  mutate(time_year = time / 365.25)
```

```{r fig2a}
sim_typ |>
  ggplot(aes(time_year, fpg, colour = cohort, linetype = cohort)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c(
    "Pioglitazone (F)" = "#1F77B4",
    "Pioglitazone (M)" = "#1F77B4",
    "Control (F)"      = "#D62728",
    "Control (M)"      = "#D62728"
  )) +
  scale_linetype_manual(values = c(
    "Pioglitazone (F)" = "solid",
    "Pioglitazone (M)" = "dashed",
    "Control (F)"      = "solid",
    "Control (M)"      = "dashed"
  )) +
  labs(x = "Time (years)",
       y = "FPG (mg/dL)",
       colour = NULL, linetype = NULL,
       title = "Figure 2a - simulated FPG trajectories",
       caption = "Replicates Figure 2a of Stringer 2015 (typical-value simulation).") +
  theme_bw() +
  theme(legend.position = "bottom")
```

```{r fig2b}
sim_typ |>
  ggplot(aes(time_year, hba1c, colour = cohort, linetype = cohort)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c(
    "Pioglitazone (F)" = "#1F77B4",
    "Pioglitazone (M)" = "#1F77B4",
    "Control (F)"      = "#D62728",
    "Control (M)"      = "#D62728"
  )) +
  scale_linetype_manual(values = c(
    "Pioglitazone (F)" = "solid",
    "Pioglitazone (M)" = "dashed",
    "Control (F)"      = "solid",
    "Control (M)"      = "dashed"
  )) +
  labs(x = "Time (years)",
       y = "HbA1c (%)",
       colour = NULL, linetype = NULL,
       title = "Figure 2b - simulated HbA1c trajectories",
       caption = "Replicates Figure 2b of Stringer 2015 (typical-value simulation).") +
  theme_bw() +
  theme(legend.position = "bottom")
```

## 3. Comparison against published summary values

Stringer 2015 reports several derived summary values in the Results
section. The table below compares the model's typical-value predictions
against those (in mixed-sex form, mirroring the paper's 60% male / 40%
female cohort composition).

```{r summary-table}
# Weight female / male typical curves per cohort using the reported 62.2%
# male composition (Table 1 pooled sex ratio 60.6-61.9% male).
pct_male <- 0.622
sim_mix <- sim_typ |>
  group_by(TRT, time) |>
  summarise(
    fpg   = sum(fpg   * ifelse(SEXF == 1, 1 - pct_male, pct_male)),
    hba1c = sum(hba1c * ifelse(SEXF == 1, 1 - pct_male, pct_male)),
    .groups = "drop"
  ) |>
  mutate(cohort = ifelse(TRT == 1, "Pioglitazone", "Control"),
         time_year = time / 365.25)

# Model summary values
extract_at <- function(df, t_days, param) {
  df |> group_by(cohort) |>
    summarise(v = approx(time, .data[[param]], xout = t_days)$y, .groups = "drop")
}

fpg_max_drop <- sim_mix |>
  group_by(cohort) |>
  summarise(delta_fpg = round(min(fpg) - fpg[time == 0], 1), .groups = "drop")

hba1c_2p5 <- extract_at(sim_mix, 2.5 * 365, "hba1c") |>
  rename(hba1c_2p5yr = v) |> mutate(hba1c_2p5yr = round(hba1c_2p5yr, 2))

fpg_5 <- extract_at(sim_mix, 5 * 365, "fpg") |>
  rename(fpg_5yr = v) |> mutate(fpg_5yr = round(fpg_5yr, 1))

model_summary <- fpg_max_drop |>
  left_join(hba1c_2p5, by = "cohort") |>
  left_join(fpg_5, by = "cohort")

published <- tibble::tribble(
  ~cohort,        ~delta_fpg_pub, ~hba1c_2p5yr_pub, ~fpg_5yr_pub,
  "Pioglitazone", -21,             NA,               147,
  "Control",       -9,             NA,               160
)
# Paper Results also report observed HbA1c at 2.5 years:
# Pioglitazone 7.3%; Control 7.8% (mean observed).
published$hba1c_2p5yr_pub[published$cohort == "Pioglitazone"] <- 7.3
published$hba1c_2p5yr_pub[published$cohort == "Control"]      <- 7.8

compare <- model_summary |>
  left_join(published, by = "cohort") |>
  dplyr::rename(
    "Treatment cohort"                = cohort,
    "Model max FPG drop (mg/dL)"      = delta_fpg,
    "Paper max FPG drop (mg/dL)"      = delta_fpg_pub,
    "Model HbA1c at 2.5 yr (%)"       = hba1c_2p5yr,
    "Paper mean HbA1c at 2.5 yr (%)"  = hba1c_2p5yr_pub,
    "Model FPG at 5 yr (mg/dL)"       = fpg_5yr,
    "Paper FPG at 5 yr (mg/dL)"       = fpg_5yr_pub
  )

knitr::kable(
  compare,
  caption = "Typical-value model predictions vs Stringer 2015 Results summary values (paper Results section).",
  align   = c("l", rep("r", 6))
)
```

The paper reports the maximum simulated FPG reductions as -21 mg/dL for
pioglitazone and -9 mg/dL for control (Results, "Drug effect model"), the
predicted median 5-year FPG values of 147 mg/dL (pioglitazone) and
~160 mg/dL (control) (Results, "Model-based simulation results"), and
mean observed HbA1c at 2.5 years of 7.3% (pioglitazone) vs 7.8% (control)
(Results, "Observed data analysis"). The typical-value simulation above
sits within a few mg/dL and a few tenths of a percent of these
values.

# Assumptions and deviations

* **Sex-mix composition for mixed-sex comparisons.** Stringer 2015 does
  not report a joint male / female typical trajectory in Figure 2, and
  Table 1 gives per-cohort male fractions of 62.8% (pioglitazone) and
  61.6% (control). We composited the male- and female-typical
  simulations at the pooled 62.2% male fraction to compare with the
  paper's median values in the Results narrative. The
  per-cohort composition differs by less than 1 percentage point, so a
  cohort-specific weighting would not materially move the compared
  values.
* **Box-Cox transformation on the BSL_HbA1c eta.** The paper describes
  the IIV on the HbA1c baseline as a Box-Cox transformation
  (lambda = 3.28) to account for right-skewness in the individual data
  (Methods, "Inter-individual variability and residual error"). nlmixr2
  does not natively encode a Box-Cox eta shape, so this implementation
  uses a log-normal eta on the baseline HbA1c as a practical
  approximation; the omega^2 value of 0.01 is preserved. This affects
  only the shape of the between-subject variability in HbA1c, not the
  typical trajectory or the covariate relations.
* **Interpretation of the OMEGA-BLOCK off-diagonal.** Table 2 reports
  "Correlation (omega^2 BSL HbA1c, omega^2 BSL FPG) = 0.01"; we
  interpreted this as the covariance element of the NONMEM OMEGA BLOCK
  (yielding a correlation coefficient of ~0.577 between the two baseline
  etas). This is the standard NONMEM reporting convention for OMEGA
  BLOCK entries; the alternative reading (0.01 as the correlation
  coefficient) would give a nearly independent structure that would
  ordinarily be reported as separate diagonal omegas rather than an
  OMEGA BLOCK.
* **KinH derivation from steady state.** KinH is not reported in
  Table 2. It is derived per subject from the HbA1c steady-state balance
  at t = 0 (with no drug and no disease progression):
  `KinH_i = (KoutH * BSL_HbA1c_i - KinZ) / BSL_FPG_i^c`. This mirrors
  the paper's cascading indirect-response construction (the analogous
  KinG is derived from the FPG steady state as `KinG = BSL_FPG * KoutG`,
  which is the form printed in Eq. 1). Sensitivity: for the reported
  typical values (KoutH = 0.072/day, BSL_HbA1c = 7.83%, KinZ = 0.28/day,
  BSL_FPG = 156 mg/dL, c = 1.91), KinH is approximately 1.84e-5
  %-HbA1c/(mg/dL)^c/day.
* **ET50 pioglitazone = 0.** With ET50 fixed to zero, the drug effect
  DEF = Emax * t / (ET50 + t) is Emax for all t > 0 (with a small
  epsilon in the denominator to avoid 0/0 at t = 0). This reproduces
  the paper's assertion that "the maximum effect of titration [was]
  achieved for drugs in the pioglitazone group by the time of first
  FPG sample collection at 3 months" (Methods, "Drug effect model").
* **Time-unit conversion for annualized disease-progression rates.**
  Stringer 2015 reports FPGDP and DPind in years^-1. The model uses
  time in days (matching the days-scale KoutG, KoutH, ET50), so both
  rates are divided by 365.25 inside `model()` to produce a consistent
  fractional increase per day.
* **No PK.** The paper explicitly states "no pharmacokinetic data were
  collected for any of the treatments" (Discussion, Limitations). The
  drug effect is therefore time-driven within cohort (Emax and ET50
  differ by treatment group), which is why the model has no dose
  compartments and no absorption / disposition parameters.
