---
title: "PD Model Selection and Comparison"
execute:
  cache: true
  error: false
engine: julia
author:
  - Vijay Ivaturi
format:
  html:
    include-in-header:
      text: |
        <style>
          .header-container {
            width: 100%;
          }
          .logo {
            display: block;
            width: 150px;
            margin-left: -2.5%;
          }
        </style>
        <link rel = "shortcut icon" href = "https://pumas.ai/favicon.ico" />
    include-before-body:
      text: |
        <div class="header-container">
          <img src="https://pumas-assets.s3.amazonaws.com/CompanyLogos/PumasAI/RGB/PNG/Pumas+AI+Primary%404x.png" alt="Pumas Logo" class="logo">
        </div>
    self-contained-math: true
    html-math-method: katex
    anchor-sections: true
    theme: default
    toc: true
    toc-depth: 4
    toc-expand: 2
    toc-location: left
    toc-title: Contents
    number-sections: true
    code-summary: Show/Hide Code
fig-format: svg
fig-width: 8
fig-height: 6
license: CC BY-SA 4.0
bibliography: references.bib
mermaid:
  theme: neutral
---

```{julia}
#| echo: false
#| output: false
using Pumas
using PumasUtilities
using PumasPlots
using CairoMakie
using AlgebraOfGraphics
using DataFrames
using DataFramesMeta
using Random
using PharmaDatasets
```

## Learning Objectives

By the end of this tutorial, you will be able to:

  - Apply a systematic approach to PD model selection
  - Use exploratory data analysis to identify model type
  - Generate and interpret hysteresis plots
  - Compare competing models statistically
  - Avoid common pitfalls in PD modeling
  - Make scientifically justified model choices

## Introduction

Selecting the appropriate PD model is crucial for accurate characterization
of drug-effect relationships [@gabrielsson2016; @danhof2007].
This tutorial provides a systematic framework
for choosing between the four major PD model paradigms:

 1. **Direct Response**: Immediate equilibrium
 2. **Effect Compartment**: Delayed equilibrium (hysteresis)
 3. **Indirect Response**: Turnover-mediated effects
 4. **K-PD**: Irreversible binding mechanisms

```{mermaid}
%%| label: fig-decision-framework
%%| fig-cap: "Decision framework for PD model selection"
flowchart TD
    Start[PK-PD Data] --> EDA[Exploratory Data Analysis]

    EDA --> Q1{Plot Effect vs.<br/>Concentration}
    Q1 --> |Single curve| Direct[Consider<br/>Direct Response]
    Q1 --> |Loop observed| Q2{Hysteresis<br/>direction?}

    Q2 --> |Counter-clockwise| EC[Consider<br/>Effect Compartment]
    Q2 --> |Clockwise| Q3{Tolerance or<br/>indirect mechanism?}

    Q3 --> |Tolerance| Tol[Tolerance Model]
    Q3 --> |Indirect| IDR[Consider IDR]

    EDA --> Q4{Effect persists<br/>after drug clears?}
    Q4 --> |Moderately| IDR
    Q4 --> |Dramatically<br/>days/weeks| KPD[Consider K-PD]

    EDA --> Q5{Known mechanism<br/>of action?}
    Q5 --> |Receptor binding| Direct
    Q5 --> |Enzyme turnover| IDR
    Q5 --> |Covalent binding| KPD

    style Direct fill:#e1f5fe
    style EC fill:#fff3e0
    style IDR fill:#e8f5e9
    style KPD fill:#fce4ec
```

## Step 1: Exploratory Data Analysis

Before fitting any model, thorough exploratory analysis is essential.

### Plotting Effect vs. Time

```{julia}
# Generate example data with different characteristics
# We'll create datasets that represent each model type

# Common PK component
pk_model = @model begin
    @param begin
        pop_Ka ∈ RealDomain(; lower = 0)
        pop_CL ∈ RealDomain(; lower = 0)
        pop_Vc ∈ RealDomain(; lower = 0)
    end
    @pre begin
        Ka = pop_Ka
        CL = pop_CL
        Vc = pop_Vc
    end
    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end
    @vars begin
        Cp = Central / Vc
    end
end

pk_params = (pop_Ka = 1.0, pop_CL = 0.5, pop_Vc = 20.0)
dosing = DosageRegimen(100, time = 0, cmt = 1)
subject = Subject(id = 1, events = dosing)
sim_pk = simobs(pk_model, subject, pk_params, obstimes = 0:0.5:24, simulate_error = false)
```

### The Hysteresis Plot

The **hysteresis plot** (effect vs. concentration) is the most informative diagnostic:

```{julia}
#| label: fig-hysteresis-types
#| fig-cap: "Different hysteresis patterns and their model implications"

# Convert simulation to DataFrame
pk_df = DataFrame(sim_pk)

# Generate three scenarios
E0, Emax, EC50 = 20.0, 80.0, 2.5
conc = pk_df.Cp
times = pk_df.time

# Direct response: single curve
effect_direct = E0 .+ (Emax .* conc) ./ (EC50 .+ conc)

# Counter-clockwise hysteresis (effect compartment)
# Simulate delayed concentration
ce = zeros(length(conc))
ke0 = 0.3
dt = 0.5
for i = 2:length(ce)
    ce[i] = ce[i-1] + dt * ke0 * (conc[i] - ce[i-1])
end
effect_ccw = E0 .+ (Emax .* ce) ./ (EC50 .+ ce)

# Clockwise hysteresis (tolerance)
# Effect decreases over time at same concentration
tolerance = 1.0 .- 0.5 .* (1 .- exp.(-0.1 .* times))
effect_cw = (E0 .+ (Emax .* conc) ./ (EC50 .+ conc)) .* tolerance

# Create DataFrames for each scenario
direct_df = DataFrame(conc = conc, effect = effect_direct, type = "direct")
ccw_df = DataFrame(conc = conc, effect = effect_ccw, type = "ccw")
cw_df = DataFrame(conc = conc, effect = effect_cw, type = "cw")

# Plot
fig = Figure(size = (900, 400))

# Direct response
direct_plt =
    data(direct_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :steelblue) +
        visual(Scatter, markersize = 8, color = :steelblue)
    )
draw!(fig[1, 1], direct_plt; axis = (title = "No Hysteresis\n(Direct Response)",))

# Counter-clockwise
ccw_plt =
    data(ccw_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :coral) +
        visual(Scatter, markersize = 8, color = :coral)
    )
ax2 = Axis(
    fig[1, 2],
    xlabel = "Concentration",
    ylabel = "Effect",
    title = "Counter-Clockwise\n(Effect Compartment)",
)
draw!(ax2, ccw_plt)
arrows!(
    ax2,
    [conc[8]],
    [effect_ccw[8]],
    [conc[12] - conc[8]],
    [effect_ccw[12] - effect_ccw[8]],
    arrowsize = 15,
    color = :red,
)

# Clockwise
cw_plt =
    data(cw_df) *
    mapping(:conc => "Concentration", :effect => "Effect") *
    (
        visual(Lines, linewidth = 2, color = :purple) +
        visual(Scatter, markersize = 8, color = :purple)
    )
ax3 = Axis(
    fig[1, 3],
    xlabel = "Concentration",
    ylabel = "Effect",
    title = "Clockwise\n(Tolerance/IDR)",
)
draw!(ax3, cw_plt)
arrows!(
    ax3,
    [conc[8]],
    [effect_cw[8]],
    [conc[4] - conc[8]],
    [effect_cw[4] - effect_cw[8]],
    arrowsize = 15,
    color = :red,
)

fig
```

### Key Diagnostic Patterns

| Observation             | Interpretation                  | Candidate Model    |
|:----------------------- |:------------------------------- |:------------------ |
| No loop (single curve)  | Immediate equilibrium           | Direct Response    |
| Counter-clockwise loop  | Effect lags concentration       | Effect Compartment |
| Clockwise loop          | Tolerance or indirect mechanism | IDR or Tolerance   |
| Effect >> drug duration | Irreversible or slow turnover   | K-PD or IDR        |

## Step 2: Mechanism-Based Selection

When pharmacological mechanism is known, use it to guide model selection
[@mould2012; @upton2014]:

:::{.panel-tabset}

### Receptor Binding

**Drug binds reversibly to receptor**

  - Rapid equilibration → **Direct Response**
  - Slow equilibration → **Effect Compartment**

Examples: Beta-blockers, opioids, benzodiazepines

### Enzyme Inhibition

**Drug inhibits enzyme reversibly**

  - Fast on/off → **Direct Response**
  - Product turnover involved → **Indirect Response**

Examples: Statins (HMG-CoA reductase), ACE inhibitors

### Signal Transduction

**Drug triggers downstream cascade**

  - Synthesis/degradation involved → **Indirect Response**
  - Second messenger turnover → **Indirect Response**

Examples: Corticosteroids, growth factors

### Covalent Binding

**Drug binds irreversibly to target**

  - Target turnover governs recovery → **K-PD**

Examples: Aspirin, PPIs, covalent kinase inhibitors

:::

## Step 3: Statistical Model Comparison

### Information Criteria

For comparing non-nested models, use information criteria:

$$
AIC = -2 \cdot \log(L) + 2k
$$

$$
BIC = -2 \cdot \log(L) + k \cdot \log(n)
$$

Where:

  - $L$ = maximum likelihood
  - $k$ = number of parameters
  - $n$ = number of observations

:::{.callout-tip title="Interpretation"}
**Lower AIC/BIC is better.** AIC favors model fit while BIC penalizes complexity more heavily.
A difference < 2 means models are essentially equivalent; 2-10 indicates moderate evidence for the better model; > 10 indicates strong evidence.
:::

### Likelihood Ratio Test

For **nested models** (one is a special case of another):

$$
LRT = -2 \cdot (\log L_{reduced} - \log L_{full})
$$

The LRT follows a chi-squared distribution with degrees of freedom equal to
the difference in parameter count.

## Case Study: Warfarin Anticoagulation

Let's work through a real example comparing models for warfarin's anticoagulant effect.

```{julia}
# Load warfarin data
warfarin_data = dataset("pumas/warfarin_pumas")

# Process data
warfarin_df = @chain warfarin_data begin
    @rtransform begin
        :id = string(:id)
        :FSZV = :wtbl / 70
        :FSZCL = (:wtbl / 70)^0.75
    end
end

# Create population for PKPD
pkpd_pop = read_pumas(
    warfarin_df;
    id = :id,
    time = :time,
    amt = :amt,
    cmt = :cmt,
    evid = :evid,
    covariates = [:sex, :wtbl, :FSZV, :FSZCL],
    observations = [:conc, :pca],
)
```

### Model 1: Direct Response

```{julia}
warf_direct = @model begin
    @metadata begin
        desc = "Warfarin Direct Response"
        timeu = u"hr"
    end

    @param begin
        pop_CL ∈ RealDomain(; lower = 0, init = 0.13)
        pop_Vc ∈ RealDomain(; lower = 0, init = 8.0)
        pop_tabs ∈ RealDomain(; lower = 0, init = 0.5)
        pop_lag ∈ RealDomain(; lower = 0, init = 0.1)
        pop_E0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_Emax ∈ RealDomain(; init = -50.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 1.0)
        pk_Ω ∈ PDiagDomain([0.04, 0.04, 0.04])
        pd_Ω ∈ PDiagDomain([0.04])
        σ_prop ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add ∈ RealDomain(; lower = 0, init = 0.1)
        σ_pca ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates FSZCL FSZV

    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
        E0 = pop_E0 * exp(pd_η[1])
        Emax = pop_Emax
        EC50 = pop_EC50
    end

    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end

    @dynamics Depots1Central1

    @vars begin
        cp = Central / Vc
        pca_pred = E0 + Emax * cp / (EC50 + cp)
    end

    @derived begin
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
        pca ~ @. Normal(pca_pred, σ_pca)
    end
end
```

### Model 2: Indirect Response (Type III)

```{julia}
warf_idr = @model begin
    @metadata begin
        desc = "Warfarin IDR Type III"
        timeu = u"hr"
    end

    @param begin
        pop_CL ∈ RealDomain(; lower = 0, init = 0.13)
        pop_Vc ∈ RealDomain(; lower = 0, init = 8.0)
        pop_tabs ∈ RealDomain(; lower = 0, init = 0.5)
        pop_lag ∈ RealDomain(; lower = 0, init = 0.1)
        pop_E0 ∈ RealDomain(; lower = 0, init = 100.0)
        pop_Emax ∈ RealDomain(; init = 1.0)
        pop_EC50 ∈ RealDomain(; lower = 0, init = 1.0)
        pop_turnover ∈ RealDomain(; lower = 0, init = 14.0)
        pk_Ω ∈ PDiagDomain([0.04, 0.04, 0.04])
        pd_Ω ∈ PDiagDomain([0.04])
        σ_prop ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add ∈ RealDomain(; lower = 0, init = 0.1)
        σ_pca ∈ RealDomain(; lower = 0, init = 5.0)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates FSZCL FSZV

    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
        E0 = pop_E0 * exp(pd_η[1])
        Emax = pop_Emax
        EC50 = pop_EC50
        kout = log(2) / pop_turnover
        kin0 = E0 * kout
    end

    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end

    @init begin
        PCA = E0
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - CL / Vc * Central
        PCA' = kin - kout * PCA
    end

    @vars begin
        # Plasma concentration
        cp = Central / Vc
        # Drug stimulation of production
        stim = Emax * cp / (EC50 + cp)
        kin = kin0 * (1 + stim)
        # PCA for plotting
        pca_pred = PCA
    end

    @derived begin
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
        pca ~ @. Normal(PCA, σ_pca)
    end
end
```

### Fit Both Models

```{julia}
#| output: false
# Fit direct response
fit_direct = fit(warf_direct, pkpd_pop, init_params(warf_direct), FOCE())

# Fit IDR
fit_idr = fit(warf_idr, pkpd_pop, init_params(warf_idr), FOCE())
```

### Compare Models

```{julia}
# AIC comparison
aic_direct = aic(fit_direct)
aic_idr = aic(fit_idr)

println("AIC Comparison:")
println("  Direct Response: $(round(aic_direct, digits=1))")
println("  IDR Model:       $(round(aic_idr, digits=1))")
println("  Difference:      $(round(aic_direct - aic_idr, digits=1))")
```

```{julia}
# BIC comparison
bic_direct = bic(fit_direct)
bic_idr = bic(fit_idr)

println("\nBIC Comparison:")
println("  Direct Response: $(round(bic_direct, digits=1))")
println("  IDR Model:       $(round(bic_idr, digits=1))")
println("  Difference:      $(round(bic_direct - bic_idr, digits=1))")
```

### Visual Comparison

```{julia}
# Inspect both models
insp_direct = inspect(fit_direct)
insp_idr = inspect(fit_idr)
```

```{julia}
#| label: fig-model-comparison-gof
#| fig-cap: "Goodness-of-fit comparison for PCA observations"
# GOF for direct model
goodness_of_fit(insp_direct; observations = [:pca])
```

```{julia}
#| label: fig-model-comparison-gof-idr
#| fig-cap: "Goodness-of-fit for IDR model"
goodness_of_fit(insp_idr; observations = [:pca])
```

## Common Pitfalls

:::{.callout-warning}

## Pitfall 1: Overparameterization

Adding parameters without data support leads to:

  - Poor parameter precision (wide confidence intervals)
  - Model instability
  - Overfitting

**Solution**: Start simple, add complexity only when data supports it.
:::

:::{.callout-warning}

## Pitfall 2: Ignoring Mechanism

Choosing models based purely on fit statistics without considering:

  - Known pharmacology
  - Biological plausibility
  - Clinical interpretation

**Solution**: Use mechanism to inform model selection, then validate with data.
:::

:::{.callout-warning}

## Pitfall 3: Sparse PD Data

With limited PD observations:

  - Hysteresis may not be detectable
  - Turnover parameters poorly estimated
  - Complex models not identifiable

**Solution**: Match model complexity to data richness.
:::

## Decision Checklist

Use this checklist when selecting a PD model:

 1. **Exploratory Analysis**
    
      + [ ] Plot effect vs. time for all subjects
      + [ ] Create hysteresis plots (effect vs. concentration)
      + [ ] Assess variability in response

 2. **Mechanistic Considerations**
    
      + [ ] Review drug's mechanism of action
      + [ ] Identify target type (receptor, enzyme, etc.)
      + [ ] Consider known pharmacology
 3. **Model Fitting**
    
      + [ ] Start with simplest appropriate model
      + [ ] Assess parameter identifiability
      + [ ] Check convergence and stability
 4. **Model Comparison**
    
      + [ ] Compare nested models with LRT
      + [ ] Compare non-nested models with AIC/BIC
      + [ ] Evaluate goodness-of-fit diagnostics
 5. **Validation**
    
      + [ ] Visual predictive checks
      + [ ] Residual analysis
      + [ ] Clinical plausibility of parameters

## Summary

In this tutorial, we covered:

 1. **Systematic approach**: EDA → Mechanism → Statistical comparison
 2. **Hysteresis analysis**: Key diagnostic for model selection
 3. **Information criteria**: AIC/BIC for model comparison
 4. **Case study**: Warfarin model comparison
 5. **Common pitfalls**: Overparameterization, ignoring mechanism, sparse data

### Key Takeaways

  - Always start with exploratory analysis
  - Let mechanism inform model choice
  - Use statistics to compare, not choose
  - Simpler models are often better
  - Parameter interpretability matters

### Quick Reference Table

| Data Pattern                  | Mechanism              | Recommended Model  |
|:----------------------------- |:---------------------- |:------------------ |
| No hysteresis                 | Rapid receptor binding | Direct Response    |
| Counter-clockwise loop        | Distribution delay     | Effect Compartment |
| Effect outlasts concentration | Turnover process       | Indirect Response  |
| Effect >> drug half-life      | Irreversible binding   | K-PD               |
| Clockwise loop                | Tolerance/adaptation   | Tolerance Model    |
