模型模拟相关函数

SimulationResult

class SimulationResult

import polars as pl
import numpy as np

from mas.model import *


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


scenario = (
    EventTable()
    .add_dosing(
        time=0,
        amt=100,
    )
    .add_sampling(time=np.arange(0, 121))
    .with_covariates(DVID=1)
    .with_ids(np.arange(0, 1001))
    .merge(
        EventTable()
        .add_sampling(
            time=np.arange(0, 121),
        )
        .with_covariates(DVID=2)
        .with_ids(np.arange(0, 1001))
    )
)

pop_model = PopModel(mod=WarfarinPkPd, data=scenario)
res = pop_model.simulate()
res

🔧 CXX compiler C:\Program Files\Maspectra\resources\static\kits\mingw64\bin\g++.EXE
📦 Compiling build target...
🔗 Linking dynamic library...
✅ Compilation Finished
[WARNING] Casting DVID to numeric dtype

Summary

N Subject	N Obs	Seed	Simulation Time
1001	243243	1234	0:00:02.57

output

property output

模拟得到的数据集。

类型：

DataFrame

示例：

res.output.head(10)

shape: (10, 28)

ID	TIME	DV	AMT	CMT	EVID	MDV	SS	RATE	ADDL	II	DVID	_IPRED_	_PRED_	_Y_	alag	ce	cl	conc	e0	ec50	eff	emax	k	ka	ke0	v	y
i64	f64	f64	f64	i64	i64	i64	i64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64
0	0.0	null	100.0	1	1	1	null	null	null	null	1.0	0.0	0.0	-0.144564	1.027434	0.0	0.131631	0.0	117.382865	9.965346	117.382865	-250.0	0.011275	2.330835	0.181635	11.674763	-0.144564
0	0.0	null	null	2	0	1	null	null	null	null	1.0	0.0	0.0	-0.201869	1.027434	0.0	0.131631	0.0	117.382865	9.965346	117.382865	-250.0	0.011275	2.330835	0.181635	11.674763	-0.201869
0	0.0	null	null	2	0	1	null	null	null	null	2.0	117.382865	100.0	122.995336	1.027434	0.0	0.131631	0.0	117.382865	9.965346	117.382865	-250.0	0.011275	2.330835	0.181635	11.674763	122.995336
0	1.0	null	null	2	0	1	null	null	null	null	1.0	0.0	2.60881	-0.317276	1.027434	0.0	0.131631	0.0	117.382865	9.965346	117.382865	-250.0	0.011275	2.330835	0.181635	11.674763	-0.317276
0	1.0	null	null	2	0	1	null	null	null	null	2.0	117.382865	98.852754	108.755962	1.027434	0.0	0.131631	0.0	117.382865	9.965346	117.382865	-250.0	0.011275	2.330835	0.181635	11.674763	108.755962
0	2.0	null	null	2	0	1	null	null	null	null	1.0	7.621264	9.864767	8.258265	1.027434	0.851663	0.131631	7.621264	117.382865	9.965346	97.69945	-250.0	0.011275	2.330835	0.181635	11.674763	8.258265
0	2.0	null	null	2	0	1	null	null	null	null	2.0	97.69945	70.495973	95.865509	1.027434	0.851663	0.131631	7.621264	117.382865	9.965346	97.69945	-250.0	0.011275	2.330835	0.181635	11.674763	95.865509
0	3.0	null	null	2	0	1	null	null	null	null	1.0	8.331094	11.578074	9.121094	1.027434	2.060493	0.131631	8.331094	117.382865	9.965346	74.548165	-250.0	0.011275	2.330835	0.181635	11.674763	9.121094
0	3.0	null	null	2	0	1	null	null	null	null	2.0	74.548165	41.10057	73.872453	1.027434	2.060493	0.131631	8.331094	117.382865	9.965346	74.548165	-250.0	0.011275	2.330835	0.181635	11.674763	73.872453
0	4.0	null	null	2	0	1	null	null	null	null	1.0	8.315001	11.869852	8.120068	1.027434	3.102918	0.131631	8.315001	117.382865	9.965346	58.023066	-250.0	0.011275	2.330835	0.181635	11.674763	8.120068

data

property data

用于模拟的原始数据集。

类型：

DataFrame

seed

property seed

设定的随机种子号。

类型：

int

elapsed_time

property elapsed_time

模拟耗时。

类型：

timedelta

plot_ind()

plot_ind(filter=None, x_name=None, y_name='_IPRED_') → Plot

绘制个体模拟数据折线图。

参数：

• filter (Expr, optional): 筛选条件，对于部分数据进行个体模拟图绘制。默认值为 None。

• x_name (str, optional)：因变量名。若为 None，则自动推断。默认值为 None。

• y_name (str, optional)：自变量名。默认值为 "_IPRED_"。

返回值：

绘制完成的个体模拟数据折线图。

返回类型：

Plot

示例：

res.plot_ind(pl.col("ID") < 10).with_facet_wrap("DVID", shared_y_axis=False)

plot_interval()

plot_interval(filter=None, x_name=None, y_name='_IPRED_') → Plot

绘制模拟数据区间图（5% 分位数至 95% 分位数）。

参数：

• filter (Expr, optional): 筛选条件，对于部分数据进行个体模拟图绘制。默认值为 None。

• x_name (str, optional)：因变量名。若为 None，则自动推断。默认值为 None。

• y_name (str, optional)：自变量名。默认值为 "_IPRED_"。

示例：

绘制在 100 mg 给药时华法林的血药浓度个体预测值 95% 置信区间图：

res.plot_interval(filter=pl.col("DVID") == 1)