案例: 华法林 PK-PD#
在 入门案例 中，我们已经完成了一次简单的群体 PK 建模。在本节中，我们将继续该教程，进行群体 PK-PD 建模，从数据探索到模型模拟的完整教学与实战。
群体 PK 建模#
在进行 PK-PD 建模之前，我们先通过群体 PK 模型来熟悉群体建模的分析流程。
加载建模库及数据处理#
首先，我们将导入 Maspectra 的建模库 mas.model。
并从 mas.datasets 中导入软件内置的华法林 PK-PD 示例数据 (warfarin)，导入 polars 以进行数据处理 (详见 polars 入门教程)。
from mas.model import *
from mas.datasets import warfarin
import polars as pl
使用 head(n=10) 函数以检查该数据的开头 10 行。
warfarin.head(n=10)
shape: (10, 11)
| ID | TIME | AMT | DV | DVID | MDV | WT | AGE | SEX | DOSE | EVID |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | f64 | i64 | i64 | f64 | i64 |
| 0 | 0.0 | 100.0 | null | 0 | 1 | 66.7 | 50 | 1 | 100.0 | 1 |
| 0 | 0.0 | null | null | 2 | 1 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 0.5 | null | 0.0 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 1.0 | null | 1.9 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 2.0 | null | 3.3 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 3.0 | null | 6.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 6.0 | null | 9.1 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 9.0 | null | 10.8 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 12.0 | null | 8.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 24.0 | null | 5.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
使用 tail(n=10) 函数检查末尾 10 行。
warfarin.tail(n=10)
shape: (10, 11)
| ID | TIME | AMT | DV | DVID | MDV | WT | AGE | SEX | DOSE | EVID |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | f64 | i64 | i64 | f64 | i64 |
| 32 | 36.0 | null | 7.7 | 1 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 36.0 | null | 27.0 | 2 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 48.0 | null | 6.9 | 1 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 48.0 | null | 24.0 | 2 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 72.0 | null | 4.4 | 1 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 72.0 | null | 23.0 | 2 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 96.0 | null | 3.5 | 1 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 96.0 | null | 20.0 | 2 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 120.0 | null | 2.5 | 1 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
| 32 | 120.0 | null | 22.0 | 2 | 0 | 62.0 | 21 | 1 | 93.0 | 0 |
数据中 DVID=0 的记录为给药记录；DVID=1 为 PK 观测记录；DVID=2 为 PD 观测记录。
我们将从原 PK-PD 数据中，筛选给药事件与 PK 观测事件，用于 PopPK 建模：
pk_data = warfarin.filter(pl.col("DVID") != 2)
pk_data.head(10)
shape: (10, 11)
| ID | TIME | AMT | DV | DVID | MDV | WT | AGE | SEX | DOSE | EVID |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | f64 | i64 | i64 | f64 | i64 |
| 0 | 0.0 | 100.0 | null | 0 | 1 | 66.7 | 50 | 1 | 100.0 | 1 |
| 0 | 0.5 | null | 0.0 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 1.0 | null | 1.9 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 2.0 | null | 3.3 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 3.0 | null | 6.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 6.0 | null | 9.1 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 9.0 | null | 10.8 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 12.0 | null | 8.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 24.0 | null | 5.6 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 36.0 | null | 4.0 | 1 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
备注
读取 CSV 文件
可使用 polars.read_csv(r'path/to/data.csv', infer_schema_length=None, null_value=".") 函数读取 CSV 数据文件，请注意添加 infer_schema_length=None 参数以设置根据全部行数推断数据列类型，null_value="." 代表将 . 作为数据缺失值的标识符。更多详情请参考［Polars 文档](https://docs.pola.rs/api/python/stable/reference/api/polars.read_csv.html)。
数据探索#
在开始建模前，首先我们将对数据进行探索 (exploratory data analysis, EDA)。
使用函数 plot_profile() 绘制药时曲线图，以检视吸收相和消除相的数据趋势，为选择基础模型提供依据。
p_profile = plot_profile(pk_data)
p_profile
设置 with_y_axis(scale="log") 可将 Y 坐标轴修改为对数刻度:
p_profile.with_y_axis(scale="log")
半对数坐标图形提示消除相对数血药浓度与时间基本呈线性关系，符合一室模型特征。
可以进一步编辑图形的样式，以保证美观:
p_profile = (
p_profile
.with_y_axis(scale="log", axis_label="Concentration (ng/mL)") # 设置 Y 轴的标题
.with_x_axis(axis_label="Time (hours)") # 设置 X 轴的标题
.with_layout(width=500, height=500) # 设置图形的宽度和高度 (500 px * 500 px)
)
p_profile
模型定义#
基于上述数据探索的结论，我们选择 mas.model 内置的血管外给药一室模型 EvOneCmtLinear 来建立华法林群体 PK 模型，步骤如下：
在 Maspectra 中定义模型，需要创建一个继承自 Module，或其他内置模型的类，以代表模型。下列代码中定义的模型命名为
WarfarinPK：class WarfarinPk(EvOneCmtLinear): ...
接下来需要在
__init__(self)中定义模型参数，包括固定效应 theta、个体间变异 omega 与个体内变异 sigma，以及这些参数的初值与边界值。我们使用
theta()定义固定效应参数，在其中使用参数bounds=(lower, upper)定义参数估计的上下边界值。一室血管外给药模型所需的基础药动学参数包括清除率（clearance，CL）、表观分布容积（apparent volume of distribution，V）、吸收速率常数（ka）与吸收延迟（alag），我们来分别定义它们：class WarfarinPk(EvOneCmtLinear): def __init__(self) -> None: self.pop_cl = theta(0.15, bounds=(0.01, 10)) self.pop_v = theta(8, bounds=(0.01, 20)) self.pop_ka = theta(0.5, bounds=(0.01, 5)) self.pop_alag = theta(0.8, bounds=(0.01, 24))
我们使用
omega()定义方差为 omega 的个体间变异参数，使用sigma()定义方差为 sigma 的个体内变异参数。我们将为每个药动学参数赋予对应的个体间变异参数，同时假设模型具有加和型与比例型的个体内变异参数。此外，如果我们需要固定某个参数值，使用参数fixed=True即可。class WarfarinPk(EvOneCmtLinear): def __init__(self) -> None: self.pop_cl = theta(0.15, bounds=(0.01, 10)) self.pop_cl = theta(0.15, bounds=(0.01, 10)) self.pop_ka = theta(0.5, bounds=(0.01, 5)) self.pop_alag = theta(0.8, bounds=(0.01, 24)) self.eta_cl = omega(0.05) self.eta_v = omega(0.08) self.eta_ka = omega(0.05) self.eta_alag = omega(0.15, fixed=True) self.eps_prop = sigma(0.01) self.eps_add = sigma(0.55)
定义模型所需的参数后，我们在 def pred(self) 内构建基于模型参数的因变量预测值模型（例如包括参数个体间变异形式、个体内变异形式和房室间的微分方程组等）。如果使用的是内置房室模型（例如本例中的
EvOneCmtLinear），可以直接调用 self.pred_micro() 或 self.pred_physio() 等方法完成对房室内药量A和血药浓度F的预测，而免去了手动书写房室间的微分方程组或动力学方程的解析解公式。例如：
class WarfarinPk(EvOneCmtLinear):
def __init__(self) -> None:
self.pop_cl = theta(0.15, bounds=(0.01, 10))
self.pop_v = theta(8, bounds=(0.01, 20))
self.pop_ka = theta(0.5, bounds=(0.01, 5))
self.pop_alag = theta(0.8, bounds=(0.01, 24))
self.eta_cl = omega(0.05)
self.eta_v = omega(0.08)
self.eta_ka = omega(0.05)
self.eta_alag = omega(0.15)
self.eps_prop = sigma(0.01)
self.eps_add = sigma(0.55)
def pred(self):
cl = self.pop_cl * exp(self.eta_cl) # 指数型个体间变异
v = self.pop_v * exp(self.eta_v)
ka = self.pop_ka * exp(self.eta_ka)
alag = self.pop_alag * exp(self.eta_alag)
pred_res = self.pred_physio(cl=cl, v=v, ka=ka, alag1=alag, s2=v) # 使用生理药动学参数预测
ipred = pred_res.F # 取出血药浓度预测值
y = ipred * (1 + self.eps_prop) + self.eps_add # 比例型与加和型个体内变异
return y
经过上述步骤，我们已经完成了对模型的定义。此时我们通过 PopModel 将模型定义与数据结合在一起，以便于在下一步中使用模型拟合数据：
model = PopModel(mod=WarfarinPk, data=pk_data)
PopModel 模型对象支持通过 to_nonmem() 函数以就 Maspectra 代码翻译为 NONMEM 控制代码:
model.to_nonmem()
; AUTHOR: Maspectra masmod program
; TIME: 2025-01-16 17:54:35
; NOTE: The translation of Mas Model to NONMEM may not be accurate, please check and run
$PROBLEM
$DATA WarfarinPk.csv IGNORE=#
$INPUT ID TIME AMT DV DVID MDV WT AGE SEX DOSE EVID
$SUBROUTINES ADVAN2 TRANS2
$PK
CL = THETA(1) * EXP(ETA(1))
V = THETA(2) * EXP(ETA(2))
KA = THETA(3) * EXP(ETA(3))
alag = THETA(4) * EXP(ETA(4))
ALAG1 = ALAG
ALAG1 = ALAG
S2 = V
$ERROR
ipred = F
y = ipred * (1 + EPS(1)) + EPS(2)
$THETA
(0.01, 0.15, 10.0); pop_cl
(0.01, 8.0, 20.0); pop_v
(0.01, 0.5, 5.0); pop_ka
(0.01, 0.8, 24.0); pop_alag
$OMEGA 0.05
$OMEGA 0.08
$OMEGA 0.05
$OMEGA 0.15
$SIGMA 0.01
$SIGMA 0.55
如果设置参数输出路径 dir, 控制文件名 control_filename，数据文件名 data_filename，可输出 NONMEM 代码及数据到对应文件夹下。
model.to_nonmem(dir=r"path/to/dir", control_filename="run1.mod", data_filename="data.csv")
可在 NONMEM 中运行该控制文件，以与 Maspectra 软件进行结果比对。
模型拟合#
使用 fit() 方法即可拟合模型，例如我们使用 FOCEi 算法进行拟合：
fit_res = model.fit(FOCEi(print=10))
➡️ Since build cache already exists, skip compile
✈️ First Order Conditional Estimation
Using BFGS...
* maxiter = 9999
* xtol = 1e-06
* ftol = 1e-06
* lower_bounds =
* upper_bounds =
* log_level = 0
* print = 10
* print_type = 0
* verbose = 1
* p_small = 0.1
* rel_step_size = 0.001
#6 ofv: 221.0960183
x: -7.2891e-02 8.3836e-02 1.0966e+00 4.6522e-02 2.1288e-01 -1.9917e-01 1.5544e+00 3.1610e-01 -1.5021e-02 -9.3995e-01
Gradients: -27.83990018 -4.529436581 -2.559537464 0.126858985 -20.11512122 -4.806058932 4.321817677 1.625650173 11.7040482 5.676701932
Funcalls: 62
#16 ofv: 218.0402149
x: -4.0068e-02 9.4113e-02 1.2621e+00 9.5160e-02 3.5770e-01 -1.4946e-01 1.4013e+00 1.5612e-01 -4.5825e-02 -9.5688e-01
Gradients: 0.1751094888 -0.04672872979 -0.04608552924 -0.05553771551 0.2142667902 0.0360657554 0.06309788614 0.006985671555 -0.1632259517 -0.1029274517
Funcalls: 156
converged = True
n_iter = 21
x_opt = -3.8206e-02 9.5404e-02 1.3195e+00 1.3128e-01 3.5607e-01 -1.5048e-01 1.4181e+00 1.1805e-01 -4.6250e-02 -9.5334e-01
f_opt = 217.9974593
fun_calls = 241
Computing Covariance...
Covariance computation succeed
查看拟合结果，输出 OFV、参数估计值、相对标准误、Shrinkage 等信息。
fit_res
Estimation Summary
| Estimation | Covariance | OFV | Cond. Number |
|---|---|---|---|
| ✅ Converged | ✅ Succeed | 217.997 | 46.818 |
| Number of Observations: | 250 | ||
| Number of Parameters: | 10 | ||
Parameter Estimation
| Parameter | Estimate | RSE(%) | Shrinkage(%) |
|---|---|---|---|
| Theta | |||
| pop_cl | 0.132 | 5.219 | |
| pop_v | 7.978 | 4.075 | |
| pop_ka | 1.354 | 31.956 | |
| pop_alag | 0.824 | 19.050 | |
| Omega(SD) | |||
| eta_cl | 0.289 | 14.185 | 0.000 |
| eta_v | 0.220 | 13.348 | 3.995 |
| eta_ka | 0.835 | 19.202 | 50.450 |
| eta_alag | 0.394 | 28.381 | 59.243 |
| Sigma(SD) | |||
| eps_prop | 0.086 | 14.238 | 14.952 |
| eps_add | 0.259 | 23.645 | 15.636 |
Information
| Number of Iterations | Time | ||
|---|---|---|---|
| Compilation | Estimation | Covariance | |
| 21 | 0:00:00.1 | 0:00:01.4 | 0:00:01.74 |
模型评价#
在模型拟合完成后，我们需要进行模型评价来确定模型拟合结果的优劣与模型稳定性。通常，使用以下方法：
拟合优度图（Goodness of Fit, GOF） 和 个体拟合图（Individual Fit）：观察模型的拟合性能，即观测值与模型预测值之间的差异大小。
可视化预测检验（Visual Predictive Check, VPC）：评估模型的预测性能和模型结构的适合性，判断模型是否能描述数据的趋势和变异。
自举法（Bootstrap）：评价模型的稳定性，检查模型是否受个体数据的影响较大。
我们首先使用 plot_gof() 绘制拟合优度图：
p_gof = fit_res.plot_gof()
p_gof
接下来，我们用 plot_predictions() 来绘制个体拟合图：
p_ind = fit_res.plot_predictions()
p_ind
通过观察以上两张诊断图形，我们发现模型群体或个体预测值与实际观测值的差异较小，可以认为模型的拟合性能较好。
我们使用 vpc() 方法可以执行可视化预测检验：
vpc_res = fit_res.vpc()
p_vpc = vpc_res.plot()
p_vpc
[VPC] Preparing data...
[VPC] Running simulation...
[VPC] Done
Running Auto-bin...
Bin search ended, result is 10
Auto-bin completed, finally use 10 bins.
观察上述图形与表格，我们发现基于当前模型参数估计值的模拟数据可以较好地描述观测数据的趋势，可以认为模型的预测性能较好、模型结构合理。
通过 bootstrap() 方法，我们可以运行自举法：
bs_res = fit_res.bootstrap(n_sample=500)
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 50
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 100
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 150
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 200
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 250
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 300
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 350
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 400
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 450
|=====|=====|=====|=====|=====|=====|=====|=====|=====|=====| 500
Bootstrap data sampling...
✅ Bootstrap complete.
Bootstrap 的结果呈现为:
bs_res
Bootstrap Summary
| Success Rate (%) | Elapsed Time |
|---|---|
| 100.000 | 0:24:57.0 |
Parameter Summary
| Parameter | Estimate | Median | 95%CI |
|---|---|---|---|
| Theta | |||
| pop_cl | 0.132 | 0.132 | 0.119, 0.147 |
| pop_v | 7.978 | 7.961 | 7.374, 8.627 |
| pop_ka | 1.354 | 1.310 | 0.688, 3.157 |
| pop_alag | 0.824 | 0.845 | 0.512, 1.275 |
| Omega(SD) | |||
| eta_cl | 0.289 | 0.279 | 0.195, 0.358 |
| eta_v | 0.220 | 0.216 | 0.153, 0.272 |
| eta_ka | 0.835 | 0.763 | 0.002, 1.311 |
| eta_alag | 0.394 | 0.361 | 0.004, 0.669 |
| Sigma(SD) | |||
| eps_prop | 0.086 | 0.085 | 0.064, 0.116 |
| eps_add | 0.259 | 0.260 | 0.007, 0.365 |
output 属性可呈现每一个重抽样数据集的拟合成功信息，以及相应的参数估计结果：
bs_res.output
shape: (500, 13)
| Bootstrap Model | Successful | OFV | pop_cl | pop_v | pop_ka | pop_alag | eta_cl | eta_v | eta_ka | eta_alag | eps_prop | eps_add |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | bool | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| 0 | true | 69.4758 | 0.131579 | 8.462114 | 4.955 | 1.452288 | 0.250414 | 0.233674 | 0.002236 | 0.019462 | 0.059671 | 0.252801 |
| 1 | true | 268.329685 | 0.138994 | 8.238129 | 1.296728 | 0.968547 | 0.268809 | 0.186015 | 0.815012 | 0.441414 | 0.081638 | 0.314727 |
| 2 | true | 169.406063 | 0.141233 | 8.213597 | 1.741106 | 0.999829 | 0.279609 | 0.219321 | 1.073348 | 0.310194 | 0.085819 | 0.246899 |
| 3 | true | 179.147362 | 0.125884 | 8.101268 | 1.475337 | 0.602399 | 0.289765 | 0.252763 | 0.489041 | 0.487545 | 0.069324 | 0.276462 |
| 4 | true | 168.537968 | 0.133721 | 8.038699 | 1.443719 | 0.66005 | 0.30242 | 0.1671 | 1.157737 | 0.46204 | 0.082276 | 0.26129 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 495 | true | 179.793889 | 0.135223 | 8.291445 | 1.445213 | 0.961974 | 0.276032 | 0.237377 | 0.802675 | 0.336041 | 0.080046 | 0.242667 |
| 496 | true | 214.70431 | 0.131074 | 7.501975 | 1.397202 | 0.758877 | 0.304111 | 0.17127 | 1.094197 | 0.364672 | 0.082188 | 0.28976 |
| 497 | true | 135.491211 | 0.11773 | 7.92955 | 0.740233 | 0.237041 | 0.208382 | 0.200522 | 0.173061 | 1.501518 | 0.080172 | 0.15368 |
| 498 | true | 174.68192 | 0.128695 | 7.840844 | 2.544863 | 1.104274 | 0.239358 | 0.201146 | 0.52756 | 0.211722 | 0.085132 | 0.211793 |
| 499 | true | 261.161659 | 0.145618 | 8.19315 | 2.757757 | 0.99981 | 0.302907 | 0.200138 | 1.232238 | 0.333624 | 0.090526 | 0.348543 |
在自举法结果中，重抽样数据集的参数估计值中位数与原参数估计值接近，且其 95% 置信区间内包括原估计值——可以认为模型的稳定性较好。
综上所述，我们认为此模型的结构合理、拟合性能较好且稳定性较强，可以进行接下来的模型模拟步骤。
模型模拟#
我们可以使用已构建的模型与参数估计结果来模拟不同场景下的血药浓度，以便于决定药物治疗剂量与制定临床个体化给药策略。在本例中，我们将基于上述模型参数估计结果模拟 100mg 与 200mg 给药剂量下华法林血药浓度群体预测值的差异。
首先，我们需要构建一个数据集用于模拟——使用模拟时间表 EventTable() 以快速创建各类模拟场景：
我们先在零时刻增加一条给药记录：
scenario_100mg = EventTable().add_dosing(id=1, time=0, amt=100)
scenario_100mg
shape: (1, 11)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | null | null | null |
我们再来增加一些观测记录，例如：0 至 48 小时，每 0.5 小时观测一次。
import numpy as np
scenario_100mg = scenario_100mg.add_sampling(id=1, time=np.arange(0, 48.1, 0.5))
scenario_100mg
shape: (98, 11)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | null | null | null |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 0.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 1.5 | null | null | null | 0 | 1 | null | null | null | null |
| … | … | … | … | … | … | … | … | … | … | … |
| 1 | 46.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 46.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 47.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 47.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 48.0 | null | null | null | 0 | 1 | null | null | null | null |
同样的，我们来创建 200mg 给药剂量的场景：
scenario_200mg = (
EventTable() # 创建模拟事件表
.add_dosing(id=1, time=0, amt=200) # 添加给药事件 (0h 给药 200mg)
.add_sampling(time=np.arange(0, 48.1, 0.5)) # 添加采样事件 (0-48h 每 0.5h 采样)
)
scenario_200mg
shape: (98, 11)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 200.0 | null | 1 | 1 | null | null | null | null |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 0.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 1.5 | null | null | null | 0 | 1 | null | null | null | null |
| … | … | … | … | … | … | … | … | … | … | … |
| 1 | 46.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 46.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 47.0 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 47.5 | null | null | null | 0 | 1 | null | null | null | null |
| 1 | 48.0 | null | null | null | 0 | 1 | null | null | null | null |
我们使用 simulate() 方法即可基于上述场景和参数估计结果来模拟血药浓度数据：
sim_res_100mg = fit_res.simulate(scenario_100mg)
sim_res_200mg = fit_res.simulate(scenario_200mg)
通过 output 属性可以获取模拟结果：
sim_res_100mg.output.head(10)
shape: (10, 20)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II | _IPRED_ | _PRED_ | _Y_ | alag | cl | ipred | ka | v | y |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| 1 | 0.0 | null | 100.0 | 1 | 1 | 1 | null | null | null | null | 0.0 | 0.0 | -0.132185 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.132185 |
| 1 | 0.0 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 0.0 | -0.083044 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.083044 |
| 1 | 0.5 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 0.0 | -0.152653 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.152653 |
| 1 | 1.0 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 2.650486 | -0.012038 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.012038 |
| 1 | 1.5 | null | null | 2 | 0 | 1 | null | null | null | null | 10.740431 | 7.466031 | 9.600547 | 1.026812 | 0.194958 | 10.740431 | 4.177999 | 7.960605 | 9.600547 |
| 1 | 2.0 | null | null | 2 | 0 | 1 | null | null | null | null | 12.121669 | 9.862466 | 12.743016 | 1.026812 | 0.194958 | 12.121669 | 4.177999 | 7.960605 | 12.743016 |
| 1 | 2.5 | null | null | 2 | 0 | 1 | null | null | null | null | 12.161339 | 11.030104 | 10.558598 | 1.026812 | 0.194958 | 12.161339 | 4.177999 | 7.960605 | 10.558598 |
| 1 | 3.0 | null | null | 2 | 0 | 1 | null | null | null | null | 12.036506 | 11.57377 | 13.592637 | 1.026812 | 0.194958 | 12.036506 | 4.177999 | 7.960605 | 13.592637 |
| 1 | 3.5 | null | null | 2 | 0 | 1 | null | null | null | null | 11.892885 | 11.800789 | 9.038686 | 1.026812 | 0.194958 | 11.892885 | 4.177999 | 7.960605 | 9.038686 |
| 1 | 4.0 | null | null | 2 | 0 | 1 | null | null | null | null | 11.748498 | 11.867316 | 12.226478 | 1.026812 | 0.194958 | 11.748498 | 4.177999 | 7.960605 | 12.226478 |
sim_res_200mg.output.head(10)
shape: (10, 20)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II | _IPRED_ | _PRED_ | _Y_ | alag | cl | ipred | ka | v | y |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| 1 | 0.0 | null | 200.0 | 1 | 1 | 1 | null | null | null | null | 0.0 | 0.0 | -0.132185 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.132185 |
| 1 | 0.0 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 0.0 | -0.083044 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.083044 |
| 1 | 0.5 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 0.0 | -0.152653 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.152653 |
| 1 | 1.0 | null | null | 2 | 0 | 1 | null | null | null | null | 0.0 | 5.300971 | -0.012038 | 1.026812 | 0.194958 | 0.0 | 4.177999 | 7.960605 | -0.012038 |
| 1 | 1.5 | null | null | 2 | 0 | 1 | null | null | null | null | 21.480863 | 14.932061 | 19.414259 | 1.026812 | 0.194958 | 21.480863 | 4.177999 | 7.960605 | 19.414259 |
| 1 | 2.0 | null | null | 2 | 0 | 1 | null | null | null | null | 24.243338 | 19.724933 | 25.404419 | 1.026812 | 0.194958 | 24.243338 | 4.177999 | 7.960605 | 25.404419 |
| 1 | 2.5 | null | null | 2 | 0 | 1 | null | null | null | null | 24.322677 | 22.060209 | 21.688026 | 1.026812 | 0.194958 | 24.322677 | 4.177999 | 7.960605 | 21.688026 |
| 1 | 3.0 | null | null | 2 | 0 | 1 | null | null | null | null | 24.073012 | 23.147541 | 27.136074 | 1.026812 | 0.194958 | 24.073012 | 4.177999 | 7.960605 | 27.136074 |
| 1 | 3.5 | null | null | 2 | 0 | 1 | null | null | null | null | 23.78577 | 23.601578 | 18.412403 | 1.026812 | 0.194958 | 23.78577 | 4.177999 | 7.960605 | 18.412403 |
| 1 | 4.0 | null | null | 2 | 0 | 1 | null | null | null | null | 23.496995 | 23.734633 | 24.435577 | 1.026812 | 0.194958 | 23.496995 | 4.177999 | 7.960605 | 24.435577 |
使用 plot_ind() 可以绘制个体拟合数据图，我们来比较两种场景下群体预测值（
_PRED_）的差异：
sim_plt_100mg = (
sim_res_100mg
.plot_ind(y_name="_PRED_") # 绘制个体预测值，y 轴为群体预测值 _PRED_
.with_y_axis(axis_label="PRED (mg/mL)") # 设置 y 轴标签
.with_x_axis(axis_label="Time (h)") # 设置 x 轴标签
)
sim_plt_100mg
sim_plt_200mg = (
sim_res_200mg.plot_ind(y_name="_PRED_").with_y_axis(axis_label="PRED (mg/mL)").with_x_axis(axis_label="Time (h)")
)
sim_plt_200mg
接下来你可以试着修改上述的模拟场景，来查看剂量或观测时间对于血药浓度的影响；也可以进一步构建更加复杂的模拟场景以熟悉 et() 的使用方法。
保存/读取结果#
我们可以通过以下代码以保存模型计算结果，请注意 dir 中指定的是保存文件夹路径:
# 保存拟合结果
fit_res.save_to(dir=r"path/to/save/directory/of/fit/")
# 读取拟合结果
fit_res = FitResult.read_from(dir=r"path/to/save/directory/of/fit/")
# 保存 Bootstrap 结果
bs_res.save_to(dir=r"path/to/save/directory/of/fit/")
# 读取 Bootstrap 结果
bs_res = BootstrapResult.read_from(dir=r"path/to/save/directory/of/fit/")
群体 PK-PD 建模#
群体 PK-PD 建模可以整合临床 PK 与 PD 的数据，串联剂量-暴露-效应的关系，定量分析药物剂量对药物最终疗效的影响，在临床试验设计和药物剂量优化中有着广泛应用。常见的 PK-PD 模型有直接效应模型、效应室模型、间接效应模型等。在本节接下来的内容中，我们将构建华法林的 PK-PD 模型。
数据探索#
与群体 PK 建模的流程一致，我们将对 PK 和 PD 的数据进行数据探索并根据数据的趋势来决定 PK-PD 模型结构。在此前 PK 建模的过程中，我们已经知晓了华法林血药浓度数据符合一室模型的基本特征，接下来我们对 PD 数据进行检视：
pd_data = warfarin.filter(pl.col("DVID") == 2)
pd_data.head(10)
shape: (10, 11)
| ID | TIME | AMT | DV | DVID | MDV | WT | AGE | SEX | DOSE | EVID |
|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | f64 | i64 | i64 | f64 | i64 |
| 0 | 0.0 | null | null | 2 | 1 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 24.0 | null | 44.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 36.0 | null | 27.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 48.0 | null | 28.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 72.0 | null | 31.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 96.0 | null | 60.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 120.0 | null | 65.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 0 | 144.0 | null | 71.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 1 | 0.0 | null | 100.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
| 1 | 24.0 | null | 49.0 | 2 | 0 | 66.7 | 50 | 1 | 100.0 | 0 |
p_profile_pd = plot_profile(pd_data)
p_profile_pd
基于上述图形，我们可以推测 PK 与 PD 间可能存在滞后效应。plot_bivariate() 可以帮助我们进一步明确同一时间点下 PK 与 PD 数据之间的关系：
p_biv = (
plot_bivariate(warfarin)
.with_x_axis(axis_label="Concentration (mg/mL)")
.with_y_axis(axis_label="PCA")
)
p_biv
观察上图，我们可以发现 PK 与 PD 数据间存在明显的滞后效应，因此我们将选择间接效应模型进行建模。
模型构建#
在效应室模型中，我们需要额外定义周转室 (turnover compartment) 并编写房室间的微分方程组。可以使用 OdeModule 和 compartment 来构建常微分方程模型和定义新的房室。
class WarfarinPkPd(OdeModule):
"""warfarin PK-PD model"""
def __init__(self):
# LSODA 对应 NONMEM 的 ADVAN13 子程序
# 该算法针对微分方程刚性问题，提供动态步长调整，在复杂模型中计算效果较好
super().__init__(solver=odeint.LSODA(rel_tol=1e-9, abs_tol=1e-9))
# 定义房室
self.cmt_depot = compartment(default_dose=True) # Depot Compartment
self.cmt_central = compartment(default_obs=True) # Central Compartment
self.cmt_turnover = compartment() # Turnover Compartment
# 定义 PK 参数
self.pop_cl = theta(0.131, bounds=(0.01, 1), fixed=True) # CL (L/h/70kg)
self.pop_v = theta(8.06, bounds=(0.01, 20), fixed=True) # V (L/70kg)
self.pop_ka = theta(1.325, bounds=(0.01, None), fixed=True) # KA (1/h)
self.pop_alag = theta(1.03, bounds=(0.01, 24), fixed=True) # ALAG1 (h)
# 定义 PD 参数
self.pop_ebsl = theta(95.228, bounds=(0.01, 200)) # EBSL
self.pop_emax = theta(-1.0, fixed=True) # EMAX
self.pop_ec50 = theta(1.59597, bounds=(0.01, 10)) # EC50
self.pop_kout = theta(0.0349133, bounds=(0.01, 100)) # KOUT
# 定义 PK 参数个体间变异 (IIV)
self.eta_cl = omega(0.0662868, fixed=True) # ETA1[CL]
self.eta_v = omega(0.016744, fixed=True) # ETA2[V]
self.eta_ka = omega(0.452064, fixed=True) # ETA3[KA]
self.eta_alag = omega(0.302627, fixed=True) # ETA4[ALAG1]
# 定义 PD 参数个体间变异 (IIV)
self.eta_ebsl = omega(0.00203239) # ETA5[EBSL]
self.eta_emax = omega(0, fixed=True) # ETA6[EMAX]
self.eta_ec50 = omega(0.153588) # ETA7[EC50]
self.eta_kout = omega(0.0336492, ) # ETA8[KOUT]
# 定义 PK 残差变异 (RV)
self.eps_pk_prop = sigma(0.00478, fixed=True) # EPS1[PROP PK]
self.eps_pk_add = sigma(0.135, fixed=True) # EPS2[ADD PK]
# 定义 PD 残差变异
self.eps_pd_add = sigma(13.9, ) # EPS3[ADD PD]
# 定义需要额外使用的数据列:
# - 标识变量 DVID
# - 协变量: WT 体重
self.wt = column('WT')
self.dvid = column('DVID')
def pred(self):
# 定义 pk 模型参数
cl = self.pop_cl * (self.wt / 70) ** 0.75 * exp(self.eta_cl)
v2 = self.pop_v * (self.wt / 70) * exp(self.eta_v)
ka = self.pop_ka * exp(self.eta_ka)
alag1 = self.pop_alag * exp(self.eta_alag)
s2 = v2
k20 = cl / v2
# 设置 depot 房室的给药吸收延迟参数
self.cmt_depot.alag = alag1
# 定义 pd emax 模型参数
ebsl = self.pop_ebsl * exp(self.eta_ebsl)
emax = self.pop_emax * exp(self.eta_emax)
ec50 = self.pop_ec50 * exp(self.eta_ec50)
self.cmt_turnover.init_value = ebsl
kout = self.pop_kout * exp(self.eta_kout)
kin = ebsl * kout
dcp = self.cmt_central.A / s2
eff = 1 + emax * dcp / (ec50 + dcp)
# 定义微分方程
self.cmt_depot.dAdt = -ka * self.cmt_depot.A
self.cmt_central.dAdt = ka * self.cmt_depot.A - k20 * self.cmt_central.A
self.cmt_turnover.dAdt = kin * eff - kout * self.cmt_turnover.A
# 计算中央室浓度
cp = self.cmt_central.A / s2
# 计算效应室浓度
pca = self.cmt_turnover.A
if (self.dvid <= 1):
# 针对 pk 指定残差模型
ipred = cp
y = cp * (1 + self.eps_pk_prop) + self.eps_pk_add
else:
# 针对 pd 指定残差模型
ipred = pca
y = pca + self.eps_pd_add
return y
在上述模型中我们可以发现从 OdeModule 构建模型的流程与从内置模型库构建 PK 模型的流程基本一致，但我们还有如下要点需要注意：
如果需要修改常微分方程求解方法或其收敛界值，可以使用
super().__init__(odeint.DVERK(tol=tolerance))。通过
compartment()可以定义一个房室。若在其中指定default_dose=True则意为此房室为默认给药室；若default_obs=True则为默认观测室。此外，通过init_value、alag和fraction等属性可以修改房室的初始药量、吸收延迟和生物利用度。对于已经定义的房室，可以使用A属性取出房室内药量；使用dAdt定义房室内药量关于时间的微分方程。使用
column("列名")方法我们可以创建一个变量来代表数据集中的某一列，这在构建协变量模型或需要根据数据值撰写判断语句时十分有用。
模型拟合#
pkpd_model = PopModel(mod=WarfarinPkPd, data=warfarin)
fit_res_pkpd = pkpd_model.fit(FOCEi(print=10), cov=False)
[WARNING] Casting DVID to numeric dtype
➡️ Since build cache already exists, skip compile
✈️ First Order Conditional Estimation
Using BFGS...
* maxiter = 9999
* xtol = 1e-06
* ftol = 1e-06
* lower_bounds =
* upper_bounds =
* log_level = 0
* print = 10
* print_type = 0
* verbose = 1
* p_small = 0.1
* rel_step_size = 0.001
#6 ofv: 1524.440893
x: 1.2807e-01 -1.1201e-01 6.6152e-01 3.0286e-01 1.9724e-01 -7.0634e-01 1.0281e-01
Gradients: 14.14091349 32.57308862 10.17788024 1.374424811 -5.459853608 -6.627805275 -17.17337527
Funcalls: 50
converged = True
n_iter = 16
x_opt = 1.2684e-01 -2.5466e-01 6.5713e-01 2.5909e-01 2.2676e-01 -4.6255e-01 1.1855e-01
f_opt = 1521.036684
fun_calls = 144
查看拟合结果:
fit_res_pkpd
Estimation Summary
| Estimation | Covariance | OFV | Cond. Number |
|---|---|---|---|
| ✅ Converged | 1521.037 | ||
| Number of Observations: | 483 | ||
| Number of Parameters: | 7 | ||
Parameter Estimation
| Parameter | Estimate | RSE(%) | Shrinkage(%) |
|---|---|---|---|
| Theta | |||
| pop_cl | 0.131 | ||
| pop_v | 8.060 | ||
| pop_ka | 1.325 | ||
| pop_alag | 1.030 | ||
| pop_ebsl | 96.568 | ||
| pop_emax | -1.000 | ||
| pop_ec50 | 1.178 | ||
| pop_kout | 0.053 | ||
| Omega(SD) | |||
| eta_cl | 0.257 | 0.000 | |
| eta_v | 0.129 | 1.725 | |
| eta_ka | 0.672 | 47.748 | |
| eta_alag | 0.550 | 64.821 | |
| eta_ebsl | 0.053 | 12.201 | |
| eta_emax | 0.000 | 0.000 | |
| eta_ec50 | 0.445 | 3.274 | |
| eta_kout | 0.105 | 20.886 | |
| Sigma(SD) | |||
| eps_pk_prop | 0.070 | 10.626 | |
| eps_pk_add | 0.332 | 0.000 | |
| eps_pd_add | 3.373 | 8.192 | |
Information
| Number of Iterations | Time | ||
|---|---|---|---|
| Compilation | Estimation | Covariance | |
| 16 | 0:00:00.1 | 0:00:49.92 | |
模型评价#
绘制 PD 数据的 GOF 图需要筛选出所有 DVID 列为 2 的数据:
p_gof_pd = fit_res_pkpd.plot_gof(filter=pl.col("DVID") == 2)
p_gof_pd
检视 PD 个体拟合图:
p_gof_pd = fit_res_pkpd.plot_predictions(filter=pl.col("DVID") == 2)
p_gof_pd
对拟合结果进行 1000 次 VPC 模拟:
vpc_pkpd = fit_res_pkpd.vpc(n_sample=1000)
[VPC] Preparing data...
[VPC] Running simulation...
[VPC] Done
呈现 PK 的 VPC 图:
p_vpc_pk = vpc_pkpd.plot(filter=pl.col("DVID") == 1)
p_vpc_pk
Running Auto-bin...
Bin search ended, result is 10
Auto-bin completed, finally use 10 bins.
PD 的 VPC 图:
p_vpc_pd = vpc_pkpd.plot(filter=pl.col("DVID") == 2)
p_vpc_pd
Running Auto-bin...
Auto-bin completed, finally use 8 bins.
小技巧
通过 filter 我们可以筛选数据供 VPC 展示。
备注
同样出于时间成本上的考虑，我们暂时不对 ODE 模型运行自举法。
模型模拟#
接下来我们使用模型模拟功能来模拟给药 100 mg 时模型个体预测值的 95% 区间。
dose_record = (
# 创建模拟事件表，对应给药记录
EventTable()
# 0h 给药 100mg
.add_dosing(time=0, amt=100)
# 给药记录 (添加数据列 DVID = 0, 体重 BW = 70 kg)
.with_covariates(DVID=0, WT=70)
# 100 个体 (ID 1-100)
.with_ids(ids=np.arange(1, 101, 1)))
pk_record = (
# 创建模拟事件表，对应 PK 采样记录
EventTable()
# 0-145h 每小时采样一次
.add_sampling(time=np.arange(0, 145, 1))
# PK 采样记录 (添加数据列 DVID = 1, 体重 BW = 70 kg)
.with_covariates(DVID=1, WT=70)
# 100 个体 (ID 1-100)
.with_ids(ids=np.arange(1, 101, 1))
)
pd_record = (
# 创建模拟事件表，对应 PD 采样记录
EventTable()
# 0-145h 每小时采样一次
.add_sampling(time=np.arange(0, 145, 1))
# PD 采样记录 (添加数据列 DVID = 2, 体重 BW = 70 kg)
.with_covariates(DVID=2, WT=70)
# 100 个体 (ID 1-100)
.with_ids(ids=np.arange(1, 101, 1))
)
scenario_pkpd = (
dose_record
# 给药事件合并 PK 采样事件
.merge(pk_record)
# 再合并 PD 采样事件
.merge(pd_record)
)
scenario_pkpd
shape: (29_100, 13)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | ADDL | II | DVID | WT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | i64 | f64 | i32 | i32 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | null | null | null | 0 | 70 |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | 1 | 70 |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | 2 | 70 |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | 1 | 70 |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | 2 | 70 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 100 | 142.0 | null | null | null | 0 | 1 | null | null | null | null | 2 | 70 |
| 100 | 143.0 | null | null | null | 0 | 1 | null | null | null | null | 1 | 70 |
| 100 | 143.0 | null | null | null | 0 | 1 | null | null | null | null | 2 | 70 |
| 100 | 144.0 | null | null | null | 0 | 1 | null | null | null | null | 1 | 70 |
| 100 | 144.0 | null | null | null | 0 | 1 | null | null | null | null | 2 | 70 |
小技巧
使用 with_facet_wrap(by) 可对图形参照指定变量分层输出。
sim_res_pkpd = fit_res_pkpd.simulate(scenario_pkpd)
p_ind_pkpd = (
sim_res_pkpd
.plot_interval()
.with_facet_wrap(by="DVID")
)
p_ind_pkpd
[WARNING] Casting WT to numeric dtype
[WARNING] Casting DVID to numeric dtype
也可以通过 filter 参数指定筛选条件，以对 pk 和 pd 数据独立绘制图形。
p_ind_pkpd_pk = (
sim_res_pkpd
.plot_interval(filter=pl.col("DVID") == 1)
.with_y_axis(axis_label="Concentration (mg/L)")
.with_x_axis(axis_label="Time (h)")
)
p_ind_pkpd_pk
p_ind_pkpd_pk = (
sim_res_pkpd
.plot_interval(filter=pl.col("DVID") == 2)
.with_y_axis(axis_label="PCA")
.with_x_axis(axis_label="Time (h)")
)
p_ind_pkpd_pk