华法林 PK-PD 建模实战#
在 入门案例 中，我们已经尝试用 Masmod 完成了一次简单的群体 PK 建模。在本节中，我们将使用 Masmod 来完成一次群体 PK-PD 建模从数据探索至模型模拟的教学与实战。
群体 PK 建模#
在进行 PK-PD 建模前，我们先通过构建一个群体 PK 模型来熟悉群体建模的分析流程。
数据探索#
导入华法林 PK-PD 示例数据：
import polars as pl
from mas.datasets import datasets
from mas.model import *
warfarin_pkpd_data = datasets.warfarin
warfarin_pkpd_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.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 |
数据中
DVID=0的记录为给药记录；DVID=1为 PK 观测记录；DVID=2为 PD 观测记录。我们将取出给药与 PK 观测记录用于 PK 建模：
pk_data = warfarin_pkpd_data.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 |
我们一般通过图形的方法来检视 PK 数据，通过观测吸收相和消除相的数据趋势以选择合适的房室模型。我们可以使用 plot_profile() 的方法来检视不同时间点下血药浓度的数据趋势：
pf_plt = plot_profile(pk_data)
pf_plt
一般地，我们会观察对数血药浓度和时间的关系来确定其符合何种房室模型。我们可以通过设置 with_y_axis(scale="log") 来将坐标轴修改为对数刻度：
pf_plt.with_y_axis(scale="log")
从上述图形我们可以发现，数据符合一室血管外给药线性消除模型的特征（消除相对数血药浓度与时间基本呈线性关系），且带有吸收延迟。
模型构建#
根据上述数据探索的结论，我们可以选择内置的一房室血管外给药模型 EvOneCmtLinear 来构建华法林群体 PK 模型。
若要在 Masmod 中构建模型，我们需要创建一个继承自 Module 或其他内置模型的类来代表模型，例如：
class WarfarinPk(EvOneCmtLinear):
pass
接下来，我们需要在
__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.theta_cl = theta(0.15, bounds=(0.01, 10))
self.theta_v = theta(8, bounds=(0.01, 20))
self.theta_ka = theta(0.5, bounds=(0.01, 5))
self.theta_alag = theta(0.8, bounds=(0.01, 24))
我们使用
omega()定义方差为 omega 的个体间变异参数，使用sigma()定义方差为 sigma 的个体内变异参数。我们将为每个药动学参数赋予对应的个体间变异参数，同时假设模型具有加和型与比例型的个体内变异参数。此外，如果我们需要固定某个参数值，使用参数fixed=True即可。
class WarfarinPk(EvOneCmtLinear):
def __init__(self) -> None:
self.theta_cl = theta(0.15, bounds=(0.01, 10))
self.theta_cl = theta(0.15, bounds=(0.01, 10))
self.theta_ka = theta(0.5, bounds=(0.01, 5))
self.theta_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.theta_cl = theta(0.15, bounds=(0.01, 10))
self.theta_v = theta(8, bounds=(0.01, 20))
self.theta_ka = theta(0.5, bounds=(0.01, 5))
self.theta_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.theta_cl * exp(self.eta_cl) # 指数型个体间变异
v = self.theta_v * exp(self.eta_v)
ka = self.theta_ka * exp(self.eta_ka)
alag = self.theta_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)
模型拟合#
使用 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
#0 ofv: 364.3540522
x: 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01 1.0000e-01
Gradients: 101.087 7.44835 -104.527 30.9969 -37.0623 16.5508 -45.9776 -13.8588 32.8152 114.414
Funcalls: 8
#1 ofv: 299.553946
x: -1.6505e-01 8.0470e-02 3.7407e-01 1.8725e-02 1.9718e-01 5.6603e-02 2.2056e-01 1.3634e-01 1.3957e-02 -2.0000e-01
Gradients: -114.076 -2.64763 -72.1449 10.9277 -34.0078 15.7933 -36.542 -12.0311 38.902 85.3421
Funcalls: 17
#11 ofv: 218.3357352
x: -4.3254e-02 9.8812e-02 1.1537e+00 1.2985e-01 3.2023e-01 -1.5320e-01 1.4055e+00 1.1531e-01 -4.2353e-02 -9.2249e-01
Gradients: -3.03182 1.55191 -1.94733 2.47344 -4.5266 -0.365343 0.687803 -0.271938 4.14002 3.55274
Funcalls: 108
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
fit_res
Estimation Summary
| OFV | Converged | Total Iterations | Compilation Time | Estimation Time | Covariance Time |
|---|---|---|---|---|---|
| 217.997 | True ✅ | 21 | 0:00:00.2 | 0:00:01.28 | 0:00:01.76 |
Parameter Estimation
| Parameter | Estimate | Shrinkage(%) | RSE(%) |
|---|---|---|---|
| Theta | |||
| theta_cl | 0.132 | 5.217 | |
| theta_v | 7.978 | 4.075 | |
| theta_ka | 1.354 | 31.821 | |
| theta_alag | 0.824 | 18.984 | |
| Omega(SD) | |||
| eta_cl | 0.289 | 0.000 | 14.191 |
| eta_v | 0.220 | 3.995 | 13.348 |
| eta_ka | 0.835 | 50.450 | 19.130 |
| eta_alag | 0.394 | 59.243 | 28.304 |
| Sigma(SD) | |||
| eps_prop | 0.086 | 14.952 | 14.138 |
| eps_add | 0.259 | 15.636 | 23.595 |
通过 plot_iter() 我们可以绘制迭代历史图，观察参数估计值是否已经收敛：
iter_plt = fit_res.plot_iterations()
iter_plt
观察上述图形，我们可以认为参数估计已经收敛，可对模型进行进一步的评价。
模型评价#
在模型拟合完成后，我们需要进行模型评价来确定模型拟合结果的优劣与模型稳定性。一般的，我们使用拟合优度图（goodness of fit, GOF）和个体拟合（individual fit）图来观察模型的拟合性能（观测值与模型预测值之间的差异大小）；使用可视化预测检验（visual predictive check, VPC）来评估模型的预测性能与模型结构是否合适（能否描述数据的趋势与其中的变异）；使用自举法（bootstrap）评价模型的稳定性（是否受个体数据的影响较大）。
我们首先使用 plot_gof() 绘制拟合优度图：
gof_plt = fit_res.plot_gof()
gof_plt
接下来，我们用 plot_ind() 来绘制个体拟合图：
ind_plt = fit_res.plot_predictions()
ind_plt
通过观察以上两张诊断图形，我们发现模型群体或个体预测值与实际观测值的差异较小，可以认为模型的拟合性能较好。
我们使用 vpc() 方法可以执行可视化预测检验：
vpc_res = fit_res.vpc()
vpc_res.plot()
[VPC] Preparing data...
[VPC] Running simulation...
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)
Bootstrap data sampling...
➡️ Bootstrap #1
➡️ Bootstrap #2
➡️ Bootstrap #3
➡️ Bootstrap #4
➡️ Bootstrap #5
➡️ Bootstrap #6
➡️ Bootstrap #7
➡️ Bootstrap #8
➡️ Bootstrap #9
➡️ Bootstrap #10
➡️ Bootstrap #11
➡️ Bootstrap #12
➡️ Bootstrap #13
➡️ Bootstrap #14
➡️ Bootstrap #15
➡️ Bootstrap #16
➡️ Bootstrap #17
➡️ Bootstrap #18
➡️ Bootstrap #19
➡️ Bootstrap #20
➡️ Bootstrap #21
➡️ Bootstrap #22
➡️ Bootstrap #23
➡️ Bootstrap #24
➡️ Bootstrap #25
➡️ Bootstrap #26
➡️ Bootstrap #27
➡️ Bootstrap #28
➡️ Bootstrap #29
➡️ Bootstrap #30
➡️ Bootstrap #31
➡️ Bootstrap #32
➡️ Bootstrap #33
➡️ Bootstrap #34
➡️ Bootstrap #35
➡️ Bootstrap #36
➡️ Bootstrap #37
➡️ Bootstrap #38
➡️ Bootstrap #39
➡️ Bootstrap #40
➡️ Bootstrap #41
➡️ Bootstrap #42
➡️ Bootstrap #43
➡️ Bootstrap #44
➡️ Bootstrap #45
➡️ Bootstrap #46
➡️ Bootstrap #47
➡️ Bootstrap #48
➡️ Bootstrap #49
➡️ Bootstrap #50
➡️ Bootstrap #51
➡️ Bootstrap #52
➡️ Bootstrap #53
➡️ Bootstrap #54
➡️ Bootstrap #55
➡️ Bootstrap #56
➡️ Bootstrap #57
➡️ Bootstrap #58
➡️ Bootstrap #59
➡️ Bootstrap #60
➡️ Bootstrap #61
➡️ Bootstrap #62
➡️ Bootstrap #63
➡️ Bootstrap #64
➡️ Bootstrap #65
➡️ Bootstrap #66
➡️ Bootstrap #67
➡️ Bootstrap #68
➡️ Bootstrap #69
➡️ Bootstrap #70
➡️ Bootstrap #71
➡️ Bootstrap #72
➡️ Bootstrap #73
➡️ Bootstrap #74
➡️ Bootstrap #75
➡️ Bootstrap #76
➡️ Bootstrap #77
➡️ Bootstrap #78
➡️ Bootstrap #79
➡️ Bootstrap #80
➡️ Bootstrap #81
➡️ Bootstrap #82
➡️ Bootstrap #83
➡️ Bootstrap #84
➡️ Bootstrap #85
➡️ Bootstrap #86
➡️ Bootstrap #87
➡️ Bootstrap #88
➡️ Bootstrap #89
➡️ Bootstrap #90
➡️ Bootstrap #91
➡️ Bootstrap #92
➡️ Bootstrap #93
➡️ Bootstrap #94
➡️ Bootstrap #95
➡️ Bootstrap #96
➡️ Bootstrap #97
➡️ Bootstrap #98
➡️ Bootstrap #99
➡️ Bootstrap #100
➡️ Bootstrap #101
➡️ Bootstrap #102
➡️ Bootstrap #103
➡️ Bootstrap #104
➡️ Bootstrap #105
➡️ Bootstrap #106
➡️ Bootstrap #107
➡️ Bootstrap #108
➡️ Bootstrap #109
➡️ Bootstrap #110
➡️ Bootstrap #111
➡️ Bootstrap #112
➡️ Bootstrap #113
➡️ Bootstrap #114
➡️ Bootstrap #115
➡️ Bootstrap #116
➡️ Bootstrap #117
➡️ Bootstrap #118
➡️ Bootstrap #119
➡️ Bootstrap #120
➡️ Bootstrap #121
➡️ Bootstrap #122
➡️ Bootstrap #123
➡️ Bootstrap #124
➡️ Bootstrap #125
➡️ Bootstrap #126
➡️ Bootstrap #127
➡️ Bootstrap #128
➡️ Bootstrap #129
➡️ Bootstrap #130
➡️ Bootstrap #131
➡️ Bootstrap #132
➡️ Bootstrap #133
➡️ Bootstrap #134
➡️ Bootstrap #135
➡️ Bootstrap #136
➡️ Bootstrap #137
➡️ Bootstrap #138
➡️ Bootstrap #139
➡️ Bootstrap #140
➡️ Bootstrap #141
➡️ Bootstrap #142
➡️ Bootstrap #143
➡️ Bootstrap #144
➡️ Bootstrap #145
➡️ Bootstrap #146
➡️ Bootstrap #147
➡️ Bootstrap #148
➡️ Bootstrap #149
➡️ Bootstrap #150
➡️ Bootstrap #151
➡️ Bootstrap #152
➡️ Bootstrap #153
➡️ Bootstrap #154
➡️ Bootstrap #155
➡️ Bootstrap #156
➡️ Bootstrap #157
➡️ Bootstrap #158
➡️ Bootstrap #159
➡️ Bootstrap #160
➡️ Bootstrap #161
➡️ Bootstrap #162
➡️ Bootstrap #163
➡️ Bootstrap #164
➡️ Bootstrap #165
➡️ Bootstrap #166
➡️ Bootstrap #167
➡️ Bootstrap #168
➡️ Bootstrap #169
➡️ Bootstrap #170
➡️ Bootstrap #171
➡️ Bootstrap #172
➡️ Bootstrap #173
➡️ Bootstrap #174
➡️ Bootstrap #175
➡️ Bootstrap #176
➡️ Bootstrap #177
➡️ Bootstrap #178
➡️ Bootstrap #179
➡️ Bootstrap #180
➡️ Bootstrap #181
➡️ Bootstrap #182
➡️ Bootstrap #183
➡️ Bootstrap #184
➡️ Bootstrap #185
➡️ Bootstrap #186
➡️ Bootstrap #187
➡️ Bootstrap #188
➡️ Bootstrap #189
➡️ Bootstrap #190
➡️ Bootstrap #191
➡️ Bootstrap #192
➡️ Bootstrap #193
➡️ Bootstrap #194
➡️ Bootstrap #195
➡️ Bootstrap #196
➡️ Bootstrap #197
➡️ Bootstrap #198
➡️ Bootstrap #199
➡️ Bootstrap #200
➡️ Bootstrap #201
➡️ Bootstrap #202
➡️ Bootstrap #203
➡️ Bootstrap #204
➡️ Bootstrap #205
➡️ Bootstrap #206
➡️ Bootstrap #207
➡️ Bootstrap #208
➡️ Bootstrap #209
➡️ Bootstrap #210
➡️ Bootstrap #211
➡️ Bootstrap #212
➡️ Bootstrap #213
➡️ Bootstrap #214
➡️ Bootstrap #215
➡️ Bootstrap #216
➡️ Bootstrap #217
➡️ Bootstrap #218
➡️ Bootstrap #219
➡️ Bootstrap #220
➡️ Bootstrap #221
➡️ Bootstrap #222
➡️ Bootstrap #223
➡️ Bootstrap #224
➡️ Bootstrap #225
➡️ Bootstrap #226
➡️ Bootstrap #227
➡️ Bootstrap #228
➡️ Bootstrap #229
➡️ Bootstrap #230
➡️ Bootstrap #231
➡️ Bootstrap #232
➡️ Bootstrap #233
➡️ Bootstrap #234
➡️ Bootstrap #235
➡️ Bootstrap #236
➡️ Bootstrap #237
➡️ Bootstrap #238
➡️ Bootstrap #239
➡️ Bootstrap #240
➡️ Bootstrap #241
➡️ Bootstrap #242
➡️ Bootstrap #243
➡️ Bootstrap #244
➡️ Bootstrap #245
➡️ Bootstrap #246
➡️ Bootstrap #247
➡️ Bootstrap #248
➡️ Bootstrap #249
➡️ Bootstrap #250
➡️ Bootstrap #251
➡️ Bootstrap #252
➡️ Bootstrap #253
➡️ Bootstrap #254
➡️ Bootstrap #255
➡️ Bootstrap #256
➡️ Bootstrap #257
➡️ Bootstrap #258
➡️ Bootstrap #259
➡️ Bootstrap #260
➡️ Bootstrap #261
➡️ Bootstrap #262
➡️ Bootstrap #263
➡️ Bootstrap #264
➡️ Bootstrap #265
➡️ Bootstrap #266
➡️ Bootstrap #267
➡️ Bootstrap #268
➡️ Bootstrap #269
➡️ Bootstrap #270
➡️ Bootstrap #271
➡️ Bootstrap #272
➡️ Bootstrap #273
➡️ Bootstrap #274
➡️ Bootstrap #275
➡️ Bootstrap #276
➡️ Bootstrap #277
➡️ Bootstrap #278
➡️ Bootstrap #279
➡️ Bootstrap #280
➡️ Bootstrap #281
➡️ Bootstrap #282
➡️ Bootstrap #283
➡️ Bootstrap #284
➡️ Bootstrap #285
➡️ Bootstrap #286
➡️ Bootstrap #287
➡️ Bootstrap #288
➡️ Bootstrap #289
➡️ Bootstrap #290
➡️ Bootstrap #291
➡️ Bootstrap #292
➡️ Bootstrap #293
➡️ Bootstrap #294
➡️ Bootstrap #295
➡️ Bootstrap #296
➡️ Bootstrap #297
➡️ Bootstrap #298
➡️ Bootstrap #299
➡️ Bootstrap #300
➡️ Bootstrap #301
➡️ Bootstrap #302
➡️ Bootstrap #303
➡️ Bootstrap #304
➡️ Bootstrap #305
➡️ Bootstrap #306
➡️ Bootstrap #307
➡️ Bootstrap #308
➡️ Bootstrap #309
➡️ Bootstrap #310
➡️ Bootstrap #311
➡️ Bootstrap #312
➡️ Bootstrap #313
➡️ Bootstrap #314
➡️ Bootstrap #315
➡️ Bootstrap #316
➡️ Bootstrap #317
➡️ Bootstrap #318
➡️ Bootstrap #319
➡️ Bootstrap #320
➡️ Bootstrap #321
➡️ Bootstrap #322
➡️ Bootstrap #323
➡️ Bootstrap #324
➡️ Bootstrap #325
➡️ Bootstrap #326
➡️ Bootstrap #327
➡️ Bootstrap #328
➡️ Bootstrap #329
➡️ Bootstrap #330
➡️ Bootstrap #331
➡️ Bootstrap #332
➡️ Bootstrap #333
➡️ Bootstrap #334
➡️ Bootstrap #335
➡️ Bootstrap #336
➡️ Bootstrap #337
➡️ Bootstrap #338
➡️ Bootstrap #339
➡️ Bootstrap #340
➡️ Bootstrap #341
➡️ Bootstrap #342
➡️ Bootstrap #343
➡️ Bootstrap #344
➡️ Bootstrap #345
➡️ Bootstrap #346
➡️ Bootstrap #347
➡️ Bootstrap #348
➡️ Bootstrap #349
➡️ Bootstrap #350
➡️ Bootstrap #351
➡️ Bootstrap #352
➡️ Bootstrap #353
➡️ Bootstrap #354
➡️ Bootstrap #355
➡️ Bootstrap #356
➡️ Bootstrap #357
➡️ Bootstrap #358
➡️ Bootstrap #359
➡️ Bootstrap #360
➡️ Bootstrap #361
➡️ Bootstrap #362
➡️ Bootstrap #363
➡️ Bootstrap #364
➡️ Bootstrap #365
➡️ Bootstrap #366
➡️ Bootstrap #367
➡️ Bootstrap #368
➡️ Bootstrap #369
➡️ Bootstrap #370
➡️ Bootstrap #371
➡️ Bootstrap #372
➡️ Bootstrap #373
➡️ Bootstrap #374
➡️ Bootstrap #375
➡️ Bootstrap #376
➡️ Bootstrap #377
➡️ Bootstrap #378
➡️ Bootstrap #379
➡️ Bootstrap #380
➡️ Bootstrap #381
➡️ Bootstrap #382
➡️ Bootstrap #383
➡️ Bootstrap #384
➡️ Bootstrap #385
➡️ Bootstrap #386
➡️ Bootstrap #387
➡️ Bootstrap #388
➡️ Bootstrap #389
➡️ Bootstrap #390
➡️ Bootstrap #391
➡️ Bootstrap #392
➡️ Bootstrap #393
➡️ Bootstrap #394
➡️ Bootstrap #395
➡️ Bootstrap #396
➡️ Bootstrap #397
➡️ Bootstrap #398
➡️ Bootstrap #399
➡️ Bootstrap #400
➡️ Bootstrap #401
➡️ Bootstrap #402
➡️ Bootstrap #403
➡️ Bootstrap #404
➡️ Bootstrap #405
➡️ Bootstrap #406
➡️ Bootstrap #407
➡️ Bootstrap #408
➡️ Bootstrap #409
➡️ Bootstrap #410
➡️ Bootstrap #411
➡️ Bootstrap #412
➡️ Bootstrap #413
➡️ Bootstrap #414
➡️ Bootstrap #415
➡️ Bootstrap #416
➡️ Bootstrap #417
➡️ Bootstrap #418
➡️ Bootstrap #419
➡️ Bootstrap #420
➡️ Bootstrap #421
➡️ Bootstrap #422
➡️ Bootstrap #423
➡️ Bootstrap #424
➡️ Bootstrap #425
➡️ Bootstrap #426
➡️ Bootstrap #427
➡️ Bootstrap #428
➡️ Bootstrap #429
➡️ Bootstrap #430
➡️ Bootstrap #431
➡️ Bootstrap #432
➡️ Bootstrap #433
➡️ Bootstrap #434
➡️ Bootstrap #435
➡️ Bootstrap #436
➡️ Bootstrap #437
➡️ Bootstrap #438
➡️ Bootstrap #439
➡️ Bootstrap #440
➡️ Bootstrap #441
➡️ Bootstrap #442
➡️ Bootstrap #443
➡️ Bootstrap #444
➡️ Bootstrap #445
➡️ Bootstrap #446
➡️ Bootstrap #447
➡️ Bootstrap #448
➡️ Bootstrap #449
➡️ Bootstrap #450
➡️ Bootstrap #451
➡️ Bootstrap #452
➡️ Bootstrap #453
➡️ Bootstrap #454
➡️ Bootstrap #455
➡️ Bootstrap #456
➡️ Bootstrap #457
➡️ Bootstrap #458
➡️ Bootstrap #459
➡️ Bootstrap #460
➡️ Bootstrap #461
➡️ Bootstrap #462
➡️ Bootstrap #463
➡️ Bootstrap #464
➡️ Bootstrap #465
➡️ Bootstrap #466
➡️ Bootstrap #467
➡️ Bootstrap #468
➡️ Bootstrap #469
➡️ Bootstrap #470
➡️ Bootstrap #471
➡️ Bootstrap #472
➡️ Bootstrap #473
➡️ Bootstrap #474
➡️ Bootstrap #475
➡️ Bootstrap #476
➡️ Bootstrap #477
➡️ Bootstrap #478
➡️ Bootstrap #479
➡️ Bootstrap #480
➡️ Bootstrap #481
➡️ Bootstrap #482
➡️ Bootstrap #483
➡️ Bootstrap #484
➡️ Bootstrap #485
➡️ Bootstrap #486
➡️ Bootstrap #487
➡️ Bootstrap #488
➡️ Bootstrap #489
➡️ Bootstrap #490
➡️ Bootstrap #491
➡️ Bootstrap #492
➡️ Bootstrap #493
➡️ Bootstrap #494
➡️ Bootstrap #495
➡️ Bootstrap #496
➡️ Bootstrap #497
➡️ Bootstrap #498
➡️ Bootstrap #499
➡️ Bootstrap #500
✅ Bootstrap complete.
bs_res
Bootstrap Summary
| Success Rate (%) | Elapsed Time |
|---|---|
| 98.000 | 0:22:03.78 |
Parameter Summary
| Parameter | Estimate | Median | 95%CI |
|---|---|---|---|
| Theta | |||
| theta_cl | 0.132 | 0.132 | 0.119, 0.148 |
| theta_v | 7.978 | 7.963 | 7.385, 8.628 |
| theta_ka | 1.354 | 1.301 | 0.699, 3.232 |
| theta_alag | 0.824 | 0.847 | 0.516, 1.275 |
| Omega(SD) | |||
| eta_cl | 0.289 | 0.279 | 0.199, 0.354 |
| eta_v | 0.220 | 0.216 | 0.158, 0.274 |
| eta_ka | 0.835 | 0.761 | 0.002, 1.281 |
| eta_alag | 0.394 | 0.362 | 0.002, 0.660 |
| Sigma(SD) | |||
| eps_prop | 0.086 | 0.085 | 0.064, 0.116 |
| eps_add | 0.259 | 0.260 | 0.001, 0.367 |
output 属性可以让我们获知对于每一个重抽样数据集的拟合成功信息与参数估计结果：
bs_res.output
shape: (500, 13)
| Bootstrap Model | Successful | OFV | theta_cl | theta_v | theta_ka | theta_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.476757 | 0.13159 | 8.462776 | 4.904765 | 1.45178 | 0.250338 | 0.233726 | 0.000494 | 0.019699 | 0.05967 | 0.252785 |
| 1 | true | 268.329685 | 0.138994 | 8.238131 | 1.296722 | 0.968545 | 0.268809 | 0.186015 | 0.815013 | 0.441415 | 0.081638 | 0.314727 |
| 2 | true | 169.415429 | 0.141252 | 8.213942 | 1.713952 | 0.999722 | 0.281448 | 0.218222 | 1.045953 | 0.312518 | 0.085809 | 0.247111 |
| 3 | true | 179.147363 | 0.125884 | 8.101209 | 1.475195 | 0.602264 | 0.289714 | 0.252747 | 0.489163 | 0.487683 | 0.069333 | 0.276466 |
| 4 | true | 168.537968 | 0.133724 | 8.03875 | 1.443796 | 0.660109 | 0.302457 | 0.167092 | 1.157487 | 0.461924 | 0.08227 | 0.261359 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 495 | true | 180.04205 | 0.134915 | 8.282561 | 1.3729 | 0.897639 | 0.276413 | 0.238503 | 0.775879 | 0.359411 | 0.080569 | 0.23854 |
| 496 | true | 214.712676 | 0.131081 | 7.502594 | 1.396143 | 0.755749 | 0.303373 | 0.169584 | 1.092751 | 0.366566 | 0.082131 | 0.290217 |
| 497 | true | 135.495917 | 0.117785 | 7.929541 | 0.739901 | 0.242934 | 0.208377 | 0.200427 | 0.173054 | 1.502029 | 0.079996 | 0.153428 |
| 498 | true | 174.669205 | 0.128694 | 7.840943 | 2.562568 | 1.110932 | 0.240236 | 0.200326 | 0.507233 | 0.212168 | 0.085127 | 0.211766 |
| 499 | true | 257.60098 | 0.145913 | 8.179922 | 3.396554 | 1.10377 | 0.300012 | 0.198947 | 1.440844 | 0.310726 | 0.084908 | 0.371665 |
在自举法结果中，重抽样数据集的参数估计值中位数与原参数估计值接近，且其 95% 置信区间内包括原估计值——可以认为模型的稳定性较好。
综上所述，我们认为此模型的结构合理、拟合性能较好且稳定性较强，可以进行接下来的模型模拟步骤。
模型模拟#
我们可以使用已构建的模型与参数估计结果来模拟不同场景下的血药浓度，以便于决定药物治疗剂量与制定临床个体化给药策略。在本例中，我们将基于上述模型参数估计结果模拟 100mg 与 200mg 给药剂量下华法林血药浓度群体预测值的差异。
首先，我们需要构建一个数据集用于模拟——使用 et() 与
+可以快速创建各类模拟场景：
我们先在零时刻增加一条给药记录：
scenario_100mg = EventTable().add_dosing(id=1, time=0, amt=100)
scenario_100mg
shape: (1, 12)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | 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, 12)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | null | null | null | null |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 0.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 1.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| … | … | … | … | … | … | … | … | … | … | … | … |
| 1 | 46.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 46.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 47.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 47.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 48.0 | null | null | null | 0 | 1 | null | null | null | null | null |
同样的，我们来创建 200mg 给药剂量的场景：
scenario_200mg = EventTable().add_dosing(id=1, time=0, amt=200).add_sampling(time=np.arange(0, 48.1, 0.5))
scenario_200mg
shape: (98, 12)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II |
|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | f64 | i64 | f64 |
| 1 | 0.0 | null | 200.0 | null | 1 | 1 | null | null | null | null | null |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 0.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 1.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| … | … | … | … | … | … | … | … | … | … | … | … |
| 1 | 46.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 46.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 47.0 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 47.5 | null | null | null | 0 | 1 | null | null | null | null | null |
| 1 | 48.0 | null | null | null | 0 | 1 | null | null | null | null | null |
我们使用 simulate() 方法即可基于上述场景和参数估计结果来模拟血药浓度数据：
sim_res_100mg = fit_res.simulate(scenario_100mg.to_data_frame())
sim_res_200mg = fit_res.simulate(scenario_200mg.to_data_frame())
通过 output 属性可以获取模拟结果：
sim_res_100mg.output.head(10)
shape: (10, 21)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II | _IPRED_ | _PRED_ | _Y_ | alag | cl | ipred | ka | v | y |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | 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 | null | 0.0 | 0.0 | -0.132185 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.132185 |
| 1 | 0.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 0.0 | -0.083044 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.083044 |
| 1 | 0.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 0.0 | -0.152653 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.152653 |
| 1 | 1.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 2.650487 | -0.012038 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.012038 |
| 1 | 1.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 10.740437 | 7.466034 | 9.600552 | 1.026812 | 0.194958 | 10.740437 | 4.178006 | 7.960605 | 9.600552 |
| 1 | 2.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 12.12167 | 9.862469 | 12.743017 | 1.026812 | 0.194958 | 12.12167 | 4.178006 | 7.960605 | 12.743017 |
| 1 | 2.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 12.161339 | 11.030107 | 10.558598 | 1.026812 | 0.194958 | 12.161339 | 4.178006 | 7.960605 | 10.558598 |
| 1 | 3.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 12.036506 | 11.573772 | 13.592637 | 1.026812 | 0.194958 | 12.036506 | 4.178006 | 7.960605 | 13.592637 |
| 1 | 3.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 11.892885 | 11.80079 | 9.038686 | 1.026812 | 0.194958 | 11.892885 | 4.178006 | 7.960605 | 9.038686 |
| 1 | 4.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 11.748498 | 11.867317 | 12.226477 | 1.026812 | 0.194958 | 11.748498 | 4.178006 | 7.960605 | 12.226477 |
sim_res_200mg.output.head(10)
shape: (10, 21)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II | _IPRED_ | _PRED_ | _Y_ | alag | cl | ipred | ka | v | y |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | 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 | null | 0.0 | 0.0 | -0.132185 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.132185 |
| 1 | 0.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 0.0 | -0.083044 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.083044 |
| 1 | 0.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 0.0 | -0.152653 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.152653 |
| 1 | 1.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 0.0 | 5.300973 | -0.012038 | 1.026812 | 0.194958 | 0.0 | 4.178006 | 7.960605 | -0.012038 |
| 1 | 1.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 21.480874 | 14.932068 | 19.414269 | 1.026812 | 0.194958 | 21.480874 | 4.178006 | 7.960605 | 19.414269 |
| 1 | 2.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 24.24334 | 19.724939 | 25.404422 | 1.026812 | 0.194958 | 24.24334 | 4.178006 | 7.960605 | 25.404422 |
| 1 | 2.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 24.322677 | 22.060213 | 21.688027 | 1.026812 | 0.194958 | 24.322677 | 4.178006 | 7.960605 | 21.688027 |
| 1 | 3.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 24.073011 | 23.147543 | 27.136073 | 1.026812 | 0.194958 | 24.073011 | 4.178006 | 7.960605 | 27.136073 |
| 1 | 3.5 | null | null | 2 | 0 | 1 | null | null | null | null | null | 23.785769 | 23.601579 | 18.412402 | 1.026812 | 0.194958 | 23.785769 | 4.178006 | 7.960605 | 18.412402 |
| 1 | 4.0 | null | null | 2 | 0 | 1 | null | null | null | null | null | 23.496995 | 23.734633 | 24.435576 | 1.026812 | 0.194958 | 23.496995 | 4.178006 | 7.960605 | 24.435576 |
使用 plot_ind() 可以绘制个体拟合数据图，我们来比较两种场景下群体预测值（
_PRED_）的差异：
sim_plt_100mg = (
sim_res_100mg.plot_ind(y_name="_PRED_").with_y_axis(axis_label="PRED (mg/mL)").with_x_axis(axis_label="Time (h)")
)
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() 的使用方法。
群体 PK-PD 建模#
群体 PK-PD 建模可以整合临床 PK 与 PD 的数据，串联剂量-暴露-效应的关系，定量分析药物剂量对药物最终疗效的影响，在临床试验设计和药物剂量优化中有着广泛应用。常见的 PK-PD 模型有直接效应模型、效应室模型、间接效应模型等。在本节接下来的内容中，我们将构建华法林的 PK-PD 模型。
数据探索#
与群体 PK 建模的流程一致，我们将对 PK 和 PD 的数据进行数据探索并根据数据的趋势来决定 PK-PD 模型结构。在此前 PK 建模的过程中，我们已经知晓了华法林血药浓度数据符合一室模型的基本特征，接下来我们对 PD 数据进行检视：
pd_data = warfarin_pkpd_data.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 |
pf_plt_pd = plot_profile(pd_data)
pf_plt_pd
基于上述图形，我们可以推测 PK 与 PD 间可能存在滞后效应。plot_bivariate() 可以帮助我们进一步明确同一时间点下 PK 与 PD 数据之间的关系：
biv_plt = plot_bivariate(warfarin_pkpd_data).with_x_axis(axis_label="Concentration").with_y_axis(axis_label="PCA")
biv_plt
观察上图，我们可以发现 PK 与 PD 数据间存在明显的滞后效应，因此我们将选择效应室模型进行建模。
参考文献
如果需要进一步了解效应室模型，可参考：
Sheiner, L. B., Stanski, D. R., Vozeh, S., Miller, R. D., & Ham, J. (1979). Simultaneous modeling of pharmacokinetics and pharmacodynamics: application to d-tubocurarine. Clinical pharmacology and therapeutics, 25(3), 358–371.
模型构建#
在效应室模型中，我们需要额外定义效应室并编写房室间的微分方程组。我们可以使用 OdeModule 和 compartment 来构建常微分方程模型和定义新的房室。
class WarfarinPkPd(OdeModule):
def __init__(self) -> None:
super().__init__(odeint.DVERK(tol=1e-3))
self.dose_cmt = compartment(default_dose=True)
self.central_cmt = compartment(default_obs=True)
self.effect_cmt = compartment()
self.theta_cl = theta(0.132, bounds=(0.01, 1))
self.theta_v = theta(7.978, bounds=(0.01, 20))
self.theta_ka = theta(1.359, bounds=(0.01, 5))
self.theta_alag = theta(0.828, bounds=(0.01, 24))
self.theta_e0 = theta(100, bounds=(0.01, 500))
self.theta_emax = theta(-250, bounds=(-500, 0))
self.theta_ec50 = theta(10, bounds=(0.01, 20))
self.theta_ke0 = theta(0.2, bounds=(0.001, 3))
self.eta_V = omega(0.08)
self.eta_Cl = omega(0.08)
self.eta_ka = omega(0.16)
self.eta_alag = omega(0.15)
self.eta_e0 = omega(0.09)
self.eta_emax = omega(0, fixed=True)
self.eta_ec50 = omega(0.4)
self.eta_ke0 = omega(0.3)
self.eps_prop_pk = sigma(0.01)
self.eps_add_pk = sigma(0.06)
self.eps_add_pd = sigma(15)
self.dvid = column("DVID")
def pred(self):
# PD 参数
cl = self.theta_cl * exp(self.eta_Cl)
v = self.theta_v * exp(self.eta_V)
ka = self.theta_ka * exp(self.eta_ka)
alag = self.theta_alag * exp(self.eta_alag)
k = cl / v
# PD 参数
e0 = self.theta_e0 * exp(self.eta_e0)
emax = self.theta_emax * exp(self.eta_emax)
ec50 = self.theta_ec50 * exp(self.eta_ec50)
ke0 = self.theta_ke0 * exp(self.eta_ke0)
# 定义微分方程组
conc = self.central_cmt.A / v
ce = self.effect_cmt.A
self.dose_cmt.alag = alag
self.dose_cmt.dAdt = -ka * self.dose_cmt.A
self.central_cmt.dAdt = ka * self.dose_cmt.A - k * self.central_cmt.A
self.effect_cmt.dAdt = ke0 * (conc - ce)
eff = e0 + (emax * ce) / (ec50 + ce)
if self.dvid == 1: # PK 预测值
y = conc * (1 + self.eps_prop_pk) + self.eps_add_pk
else: # PD 预测值
y = eff + self.eps_add_pd
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_pkpd_data)
fit_res_pkpd = pkpd_model.fit(FOCEi(print=10), cov=False)
➡️ Since build cache already exists, skip compile
[WARNING] Casting DVID to numeric dtype
✈️ 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
#1 ofv: 1611.331432
x: 9.4672e-02 1.1253e-01 9.0893e-02 1.0358e-01 8.2533e-02 1.1621e-01 1.0513e-01 -1.5478e-01 8.9808e-02 1.0081e-01 1.1537e-01 9.9547e-02 6.6929e-02 7.2822e-02 7.0000e-01 8.6872e-02 9.9782e-02 1.1493e-01
Gradients: -0.478262 -3.5792 12.1202 1.28168 21.7904 -15.1649 -4.32787 131.451 21.3932 -0.877512 -26.9196 0.955272 58.8163 48.3891 -241.082 36.1201 4.78718 22.5492
Funcalls: 27
#11 ofv: 1374.307225
x: 1.1082e-01 -2.3188e-02 -1.9254e-01 1.4262e-02 1.7565e-02 2.3598e-01 1.4541e-01 -2.7838e+00 -3.7136e-01 1.3017e-01 1.1669e+00 1.5089e-01 -1.6922e+00 -1.1241e+00 5.8892e-01 -1.0172e-01 2.2038e-01 4.0031e-02
Gradients: 8.45159 -75.9637 -4.87083 -4.44745 -381.52 55.3295 26.6224 1.35704 -32.4731 3.45339 14.4098 -1.80983 -10.2605 -34.4111 48.7432 -20.7295 -1.69051 -48.5522
Funcalls: 172
#21 ofv: 1331.106592
x: 6.0731e-02 1.3513e-01 -5.4882e-01 -4.3754e-02 6.1534e-02 1.4221e-01 1.3480e-01 -2.6489e+00 -5.2563e-02 1.1267e-01 6.7992e-01 5.9408e-01 -1.7237e+00 -7.4325e-01 1.1947e-01 -3.9544e-02 1.2592e-01 1.1421e-01
Gradients: -21.2804 14.6297 -11.9121 4.4156 50.1271 -5.79264 -1.66586 -13.4674 9.06071 0.10892 -2.85311 4.11083 -0.252057 10.3556 35.8971 -3.00237 -2.44403 18.0064
Funcalls: 314
#31 ofv: 1304.78525
x: 6.6455e-02 1.6445e-01 5.6453e-02 1.6705e-01 5.8229e-02 6.2930e-02 2.8975e-01 -2.4428e+00 -1.4012e-01 1.0790e-01 8.2439e-01 4.4052e-02 -1.6748e+00 -9.2805e-01 -6.1076e-01 -4.5500e-02 1.8659e-01 8.5743e-02
Gradients: -18.4296 25.4401 -0.501711 4.61012 0.899624 -0.516671 -0.366673 -0.271938 0.0377055 0.124685 -0.266563 -0.177376 0.00644878 0.119092 0.0769994 0.0462626 -0.172505 -0.0730409
Funcalls: 511
#41 ofv: 1303.457083
x: 1.0203e-01 9.9450e-02 5.4497e-02 9.9529e-02 5.8108e-02 5.2838e-02 3.0616e-01 -2.4488e+00 -1.5923e-01 1.0102e-01 8.2074e-01 1.0310e-01 -1.6735e+00 -9.2789e-01 -6.1279e-01 -4.7117e-02 1.8726e-01 8.5449e-02
Gradients: 0.0629703 -0.183772 -0.0346461 0.043693 0.118794 -0.0343037 -0.0200198 -0.00926661 -0.0909827 0.0141391 0.0594255 -0.0102964 0.00634902 0.00896709 -0.01621 -0.00861394 -0.0162602 0.0458895
Funcalls: 767
converged = True
n_iter = 44
x_opt = 1.0167e-01 1.0020e-01 5.4669e-02 9.8596e-02 5.8092e-02 5.2479e-02 3.0667e-01 -2.4487e+00 -1.5763e-01 1.0079e-01 8.1875e-01 1.0470e-01 -1.6737e+00 -9.2791e-01 -6.1276e-01 -4.7170e-02 1.8768e-01 8.5287e-02
f_opt = 1303.456595
fun_calls = 930
message: Minimization success
fit_res_pkpd
Estimation Summary
| OFV | Converged | Total Iterations | Compilation Time | Estimation Time | Covariance Time |
|---|---|---|---|---|---|
| 1303.457 | True ✅ | 44 | 0:00:00.1 | 0:02:06.33 |
Parameter Estimation
| Parameter | Estimate | Shrinkage(%) | RSE(%) |
|---|---|---|---|
| Theta | |||
| theta_cl | 0.132 | ||
| theta_v | 7.979 | ||
| theta_ka | 1.315 | ||
| theta_alag | 0.827 | ||
| theta_e0 | 96.690 | ||
| theta_emax | -255.939 | ||
| theta_ec50 | 11.029 | ||
| theta_ke0 | 0.018 | ||
| Omega(SD) | |||
| eta_V | 0.219 | 3.718 | |
| eta_Cl | 0.283 | 0.000 | |
| eta_ka | 0.821 | 50.077 | |
| eta_alag | 0.389 | 58.560 | |
| eta_e0 | 0.051 | 18.762 | |
| eta_emax | 0.000 | 0.000 | |
| eta_ec50 | 0.226 | 6.799 | |
| eta_ke0 | 0.269 | 10.156 | |
| Sigma(SD) | |||
| eps_prop_pk | 0.086 | 14.095 | |
| eps_add_pk | 0.267 | 14.785 | |
| eps_add_pd | 3.816 | 18.165 | |
模型评价#
gof_plt_pkpd = fit_res_pkpd.plot_gof()
gof_plt_pkpd
vpc_pkpd = fit_res_pkpd.vpc()
[VPC] Preparing data...
[VPC] Running simulation...
vpc_pk = vpc_pkpd.plot(filter=pl.col("DVID") == 1)
vpc_pk
Running Auto-bin...
Bin search ended, result is 10
Auto-bin completed, finally use 10 bins.
vpc_pd = vpc_pkpd.plot(filter=pl.col("DVID") == 2)
vpc_pd
Running Auto-bin...
Auto-bin completed, finally use 8 bins.
小技巧
通过 filter 我们可以筛选数据供 VPC 展示。
备注
同样出于时间成本上的考虑，我们暂时不对 ODE 模型运行自举法。
模型模拟#
接下来我们使用模型模拟功能来模拟给药 100 mg 时模型个体预测值的 95% 区间。
dose_record = EventTable().add_dosing(time=0, amt=100).with_covariates(DVID=0).with_ids(ids=np.arange(1, 101, 1))
pk_record = (
EventTable().add_sampling(time=np.arange(0, 145, 1)).with_covariates(DVID=1).with_ids(ids=np.arange(1, 101, 1))
)
pd_record = (
EventTable().add_sampling(time=np.arange(0, 145, 1)).with_covariates(DVID=2).with_ids(ids=np.arange(1, 101, 1))
)
scenario_pkpd = dose_record.merge(pk_record).merge(pd_record)
scenario_pkpd
shape: (29_100, 13)
| ID | TIME | DV | AMT | CMT | EVID | MDV | SS | RATE | DUR | ADDL | II | DVID |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | f64 | f64 | f64 | i64 | i64 | i64 | i64 | f64 | f64 | i64 | f64 | i32 |
| 1 | 0.0 | null | 100.0 | null | 1 | 1 | null | null | null | null | null | 0 |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | null | 1 |
| 1 | 0.0 | null | null | null | 0 | 1 | null | null | null | null | null | 2 |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | null | 1 |
| 1 | 1.0 | null | null | null | 0 | 1 | null | null | null | null | null | 2 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 100 | 142.0 | null | null | null | 0 | 1 | null | null | null | null | null | 2 |
| 100 | 143.0 | null | null | null | 0 | 1 | null | null | null | null | null | 1 |
| 100 | 143.0 | null | null | null | 0 | 1 | null | null | null | null | null | 2 |
| 100 | 144.0 | null | null | null | 0 | 1 | null | null | null | null | null | 1 |
| 100 | 144.0 | null | null | null | 0 | 1 | null | null | null | null | null | 2 |
sim_pkpd = fit_res_pkpd.simulate(scenario_pkpd.to_data_frame())
intv_plt = sim_pkpd.plot_interval().with_facet_wrap("DVID")
[WARNING] Casting DVID to numeric dtype
小技巧
对图形对象使用 filter 则可以对分组后的数据分别绘制图片。
intv_plt_pk = sim_pkpd.plot_interval(filter=pl.col("DVID") == 1).with_y_axis(axis_label="PK IPRED")
intv_plt_pk
intv_plt_pk = sim_pkpd.plot_interval(filter=pl.col("DVID") == 2).with_y_axis(axis_label="PK IPRED")
intv_plt_pk