5. 模型拟合相关函数#

5.1. mas.model.PopModel#

class mas.model.PopModel(mod, data=None)#

群体模型。对其可以进行模型拟合与模拟等操作。

此处的“群体模型”一般指代的是用于拟合群体数据的非线性混合效应模型(继承自 Module 的模型)。

参数:

  • mod (AnyModule):模型。

  • data (polars.DataFrame | pandas.DataFrame | EventTable, optional):需要被模型拟合的数据集。默认值为 None

5.1.1. data#

property data#

获取将被模型拟合的数据集。

类型:

DataFrame

5.1.2. inits#

property inits#

获取模型参数及其初值。

类型:

PopParams

5.1.3. fit()#

fit(estimation=FOCEi, cov=True) FitResult#

拟合模型。

参数:

  • estimation (AnyEstimationSettings | Sequence[AnyEstimationSettings], optional):参数估计的算法及其相关设置。当前支持的算法包括 FOCEi (first-order conditional estimation with interaction)与 FO(first-order estimation)。若传入一个序列(Sequence)则将根据传入顺序依次以对应的算法估计参数。默认值为 FOCEi

  • cov (bool, optional):是否计算协方差矩阵。默认值为 True

返回值:

包含参数估计值的模型拟合结果。

返回类型:

FitResult

示例:

使用 FOCEi 算法拟合华法林 PK 模型:

from mas.model import *
from mas.datasets import warfarin
import polars as pl


class WarfarinPk(EvOneCmtLinear.Physio):
    def __init__(self):
        super().__init__()
        self.tv_cl = theta(0.134, bounds=(0, None))
        self.tv_v = theta(7.81, bounds=(0, None))
        self.tv_ka = theta(0.571, bounds=(0, None))
        self.tv_alag = theta(0.823, bounds=(0, None))
        self.eta_cl = omega(0.049)
        self.eta_v = omega(0.0802)
        self.eta_ka = omega(0.047)
        self.eta_alag = omega(0.156)
        self.eps_prop = sigma(0.0104)
        self.eps_add = sigma(0.554)

    def pred(self):
        cl = self.tv_cl * exp(self.eta_cl)
        v = self.tv_v * exp(self.eta_v)
        ka = self.tv_ka * exp(self.eta_ka)
        alag = self.tv_alag * exp(self.eta_alag)
        F = self.solve(cl=cl, v=v, ka=ka, alag1=alag, s2=v)
        ipred = F
        y = ipred * (1 + self.eps_prop) + self.eps_add
        return y


data = warfarin.filter(pl.col("DVID") != 2)
model = PopModel(data=data, mod=WarfarinPk)
fit_res = model.fit(FOCEi(print=10))
[MTran] Start interpreting module: WarfarinPk
[MTran] Interpretation finished in 0.064722 seconds
[MTran] Vaporization started for WarfarinPk(EvOneCmtLinear.Physio)
[MTran] Vaporization finished in 0.211242 seconds
🔧 CXX compiler C:\ProgramData\mingw64\bin\g++.EXE
📦 Compiling build target...
🔗 Linking dynamic library...
✅ Compilation Finished
✈️ 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: 347.7917575
Params:       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
NParams:      1.3400e-01  7.8100e+00  5.7100e-01  8.2300e-01  2.2136e-01  2.8320e-01  2.1679e-01  3.9497e-01  1.0198e-01  7.4431e-01
Gradients:   -1.0509e+01 -8.3374e+00 -9.6925e+01  3.0800e+01 -2.7513e+01  1.7544e+01 -3.3121e+01 -1.1511e+01  4.3184e+01  1.2218e+02
Funcalls:   8
# 1    ofv: 252.9367071
Params:       1.5161e-01  1.4094e-01  5.7596e-01 -5.1245e-02  2.3511e-01  1.3849e-02  2.6265e-01  1.5653e-01 -1.1206e-01 -5.0000e-01
NParams:      1.4110e-01  8.1364e+00  9.1906e-01  7.0748e-01  2.5338e-01  2.5982e-01  2.5509e-01  4.1794e-01  8.2493e-02  4.0849e-01
Gradients:    5.8501e+01  1.9482e+01 -1.2655e+01  4.3841e-01 -2.0968e+01  1.4030e+01 -3.3114e+01 -4.1979e+00  7.4639e+00  4.2037e+01
Funcalls:   17
# 11   ofv: 219.1455277
Params:       6.6128e-02  1.1409e-01  1.1174e+00  1.3436e-01  3.6191e-01 -1.8678e-01  1.6618e+00 -3.4404e-02 -4.8287e-02 -1.0445e+00
NParams:      1.2954e-01  7.9209e+00  1.5794e+00  8.5177e-01  2.8764e-01  2.1259e-01  1.0336e+00  3.4530e-01  8.7926e-02  2.3697e-01
Gradients:   -1.2234e+01 -8.0065e+00  1.1216e+00  2.0787e+00 -3.1722e-01 -3.4493e+00  4.2760e+00 -2.3902e+00 -2.9918e+00 -5.9325e+00
Funcalls:   109
# 21   ofv: 218.0048633
Params:       8.5523e-02  1.2083e-01  9.4660e-01  8.7367e-02  3.6555e-01 -1.5218e-01  1.4398e+00  1.1198e-01 -6.5216e-02 -9.5924e-01
NParams:      1.3207e-01  7.9744e+00  1.3314e+00  8.1267e-01  2.8868e-01  2.2007e-01  8.2776e-01  3.9973e-01  8.6450e-02  2.5807e-01
Gradients:   -4.2353e-01 -4.2524e-01 -2.0438e-01 -7.2524e-01 -6.3465e-02 -7.3090e-02 -5.6033e-02 -2.5406e-02  1.8580e-02 -1.1706e-01
Funcalls:   223
converged = True
n_iter = 25
x_opt =   8.6104e-02  1.2128e-01  9.6344e-01  1.0153e-01  3.6619e-01 -1.5174e-01  1.4490e+00  9.8417e-02 -6.5853e-02 -9.5697e-01
f_opt = 217.9974593
fun_calls = 281

Estimation Summary                                                              

────────────────────────────────────────────────────────────────────────────────
    Step         OFV     
────────────────────────────────────────────────────────────────────────────────
✅ Converged 217.997     
────────────────────────────────────────────────────────────────────────────────


Parameter Estimation                                                            

────────────────────────────────────────────────────────────────────────────────
 Parameter     Estimate   Shrinkage(%)
────────────────────────────────────────────────────────────────────────────────
Theta                                 
  tv_cl        0.132                  
  tv_v         7.978                  
  tv_ka        1.354                  
  tv_alag      0.824                  
Omega(SD)                             
  eta_cl       0.289        0.000     
  eta_v        0.220        3.994     
  eta_ka       0.835        50.450    
  eta_alag     0.394        59.243    
Sigma(SD)                             
  eps_prop     0.086        14.953    
  eps_add      0.259        15.636    
────────────────────────────────────────────────────────────────────────────────


Information                                                                     

────────────────────────────────────────────────────────────────────────────────
 Number of       Time    
 Iterations              
────────────────────────────────────────────────────────────────────────────────
25           0:00:01.73  
────────────────────────────────────────────────────────────────────────────────

Estimation finished, now start covariance calculation.
Computing Covariance...
Covariance computation succeed
vpc_res = fit_res.vpc()
vpc_res.plot()
[VPC] Preparing data...
[VPC] Running simulation...
[VPC] Done in 3.030793s
Running Auto-bin...
Bin search ended, result is 10
Auto-bin completed, finally use 10 bins.

支持的算法设置#

设置

功能

print: int = 1

print 次迭代打印一次迭代历史详情。

verbose: Literal["silent", "info", "debug", "trace"] = "info"

输出详尽等级。若为 "info" 则正常打印迭代历史与其他输出信息。若为 "silent" 则不进行任何打印。

minimize_method: Literal["nm", "powell", "bobyqa", "lbfgs", "newuoa", "bfgs"] = "lbfgs"

最小化方法。

maxiter: int = 9999

最大迭代次数。

xtol: float = 1e-6

相对参数收敛界值。

ftol: float = 1e-6

相对目标函数收敛界值。

cov_method: Literal["r", "s", "rs"] = "rs"

协方差估算方法。

map_minimizer: Literal["nm", "powell", "bobyqa", "lbfgs", "newuoa", "bfgs"] = "lbfgs"

eta 参数最小化方法(仅适用于 FOCEi 算法)。

cov_minimizer: Literal["nm", "powell", "bobyqa", "lbfgs", "newuoa", "bfgs"] = "lbfgs"

协方差最小化方法(仅适用于 FOCEi 算法)。

extra: dict[str, str | int | float | bool] = dict()

各最小化方法额外设置选项

支持的各最小化方法额外设置选项#

lbfgs#

设置

功能

lbfgs_rank: int = 20

用于方向计算(direction calculation)的迭代次数。

use_approx_eigen_scaling: bool = False

使用近似特征值(approximate eigenvalue scaling)。

ls_interp_type: int = 2

用于近似目标函数的多项式的次数。

min_ls_step_size: float = 1e-9

最小线搜索(line search)步长。

gtol: float = 1e-10

相对梯度收敛界值。

num_method: int = 1

数值微分(numerical differentiation)方法。 若为 0 则使用中心(central)差分;若为 1 则使用前向(forward)差分;若为 2 则使用 Ridders 方法。

rel_step: float = 1e-3

数值微分方法相对步长(relative step size)。

BObyQA 与 NEWUOA#

设置

功能

initial_trust_region_radius: float = 0.01

初始信赖域(trust region)大小。

final_trust_region_radius: float = 1e-6

最终信赖域大小。

Nelder-Mead#

设置

功能

alpha: float = 1.0

反射(reflection)系数。

beta: float = 0.5

扩展(contraction)系数

gamma: float = 2.0

收缩(expansion)参数。

delta: float = 0.5

回退(shrinkage)系数。

adaptive: bool = True

是否使用自适应(adaptive)参数。

5.1.4. simulate()#

simulate(data=None, param=None, seed=1234) SimulationResult#

运行模型模拟。

参数:

  • data (DataFrame | EventTable, optional):用于模拟的数据集。若为 None 则使用模型将要拟合的数据集(即在创建 PopModel 时指定的 data)。支持使用 EventTable。默认值为 None

  • params (AnyPopParams, optional):用于模拟的参数值。若为 None 则使用模型中的参数初值。默认值为 None

  • seed (int, optional):随机种子。默认值为 1234

返回值:

包含模拟数据集的模拟结果。

返回类型:

SimulationResult

示例:

模拟华法林的 PK 场景

sim_res = model.simulate(data)
sim_res.output
shape: (282, 23)
IDTIMEAMTDVDVIDMDVWTAGESEXDOSEEVIDCMTSSF_IPRED__PRED__Y_alagclipredkavy
i64f64f64f64i64i64f64i64i64f64i64i64i64f64f64f64f64f64f64f64f64f64f64
00.0100.0null0166.7501100.011null0.00.00.0-0.5504621.0256040.1805160.00.7649267.788145-0.550462
00.5null0.01066.7501100.002null0.00.00.00.5856591.0256040.1805160.00.7649267.7881450.585659
01.0null1.91066.7501100.002null0.00.01.228928-1.0242611.0256040.1805160.00.7649267.788145-1.024261
02.0null3.31066.7501100.002null6.6616216.6616216.1958577.3836951.0256040.1805166.6616210.7649267.7881457.383695
03.0null6.61066.7501100.002null9.7246069.7246068.90838611.0684451.0256040.1805169.7246060.7649267.78814511.068445
3236.0null7.71062.021193.002null5.3044645.3044646.7137285.515110.8260660.1375395.3044640.48559111.9976915.51511
3248.0null6.91062.021193.002null4.6227194.6227195.4644555.1131490.8260660.1375394.6227190.48559111.9976915.113149
3272.0null4.41062.021193.002null3.5108273.5108273.6200393.867270.8260660.1375393.5108270.48559111.9976913.86727
3296.0null3.51062.021193.002null2.6663752.6663752.3981681.8023840.8260660.1375392.6663750.48559111.9976911.802384
32120.0null2.51062.021193.002null2.0250382.0250381.5887142.9416490.8260660.1375392.0250380.48559111.9976912.941649

5.1.5. to_nonmem()#

to_nonmem(dir, control_filename=None, data_filename=None) None#

将模型翻译为 NONMEM 控制文件代码。

参数:

  • dir (PathLike):翻译结果输出路径

  • control_filename (str, optional):输出的控制文件文件名。若为 None 则采用 模型类名.mod 的命名形式。默认值为 None

  • data_filename (str, optional):输出的数据集文件文件名。若为 None 则采用 模型类名.csv 的命名形式。默认值为 None

5.2. mas.model.EstimationResult#

class mas.model.EstimationResult#

模型参数估计结果。

5.2.1. ofv#

property ofv#

当前参数估计结果对应的目标函数值。

类型:

float

5.2.2. estimates#

property estimates#

模型参数估计结果。

类型:

PopParams

5.2.3. converged#

property converged#

目标函数最小化过程中参数估计是否收敛。

类型:

bool

5.2.4. n_iter#

property n_iter#

迭代次数。

类型:

int

5.2.5. est_time#

property est_time#

参数估计耗时。

类型:

timedelta

5.3. mas.model.FitResult#

class mas.model.FitResult#

模型拟合结果。

5.3.1. ofv#

property ofv

所有最小化过程结束后参数估计结果对应的目标函数值。

类型:

float

5.3.2. estimates#

property estimates

所有最小化过程结束后得到的参数估计值。

类型:

PopParams

5.3.3. converged#

property converged

所有最小化过程结束后参数估计是否收敛。

类型:

bool

5.3.4. n_iter#

property n_iter

所有最小化过程的总迭代次数。

类型:

int

5.3.5. inits#

property inits

模型参数初值。

类型:

PopParams

5.3.6. est_sequence#

property est_sequence

各最小化过程的参数估计结果。

类型:

list[EstimationResult]

5.3.7. iterations#

property iterations

各次迭代的参数估计结果。

类型:

list[IterationResults]

5.3.8. compile_time#

property compile_time

模型编译耗时。

类型:

timedelta

5.3.9. est_time#

property est_time

参数估计耗时。

类型:

timedelta

5.3.10. cov_time#

property cov_time

协方差估计耗时。

类型:

timedelta | None

5.3.11. elapsed_time#

property elapsed_time

模型拟合所花费的总耗时。

类型:

timedelta

5.3.12. output#

property output

包含观测值、群体预测值、个体预测值与个体参数值(以最后一次最小化过程的结果为准)的数据集。

类型:

DataFrame

5.3.13. covariance#

property covariance

协方差估计结果。

类型:

CovResult

5.3.14. summary()#

summary() None

打印模型拟合结果摘要。

5.3.15. compute()#

compute(residuals=None) DataFrame

计算群体预测值(population predictions)、个体预测值(individual predictions)与残差(residuals)。

参数:

  • residuals (list[Literal["CWRES", "WRES", "CWRESI", "WRESI", "CIWRES", "IWRES", "CIWRESI", "IWRESI"]], optional):需要计算的残差类型。如果为 None 则跳过残差计算。默认值为 None

返回值:

包含指定计算内容的 DataFrame

返回类型:

DataFrame

5.3.16. ebe()#

ebe() DataFrame

计算 eta 参数的经验贝叶斯估计。

返回值:

包含各个个体经验贝叶斯估计值的 DataFrame

返回类型:

DataFrame

5.3.17. eta_shr()#

eta_shr(method='sd', ebv=False) DataFrame

计算 eta 参数的收缩量。

eta 参数收缩量公式

若为标准差形式,计算公式如下:

\[ shr = [1 - \frac{SD(eta(j))}{\sqrt{omega(j, j)}}] \times 100\% \]

若为方差形式,计算公式如下:

\[ shr = [1 - \frac{Var(eta(j))}{\sqrt{omega(j, j)}}] \times 100\% \]

其中,eta(j) 为第 j 个 eta 参数的个体经验贝叶斯估计值。

参数:

  • method (Literal["sd", "var"], optional):计算收缩量的方法。使用标准差形式("sd")或方差形式("var")。默认值为 "sd"。计算公式可见:eta 参数收缩量公式

  • ebv (bool, optional):是否计算平均经验贝叶斯方差。默认值为 False

返回值:

eta 参数收缩量的数据框。

返回类型:

DataFrame

示例:

fit_res.eta_shr("sd")
shape: (1, 4)
eta_cleta_veta_kaeta_alag
f64f64f64f64
1.0000e-100.0399360.5045040.592428

5.3.18. eps_shr()#

eps_shr(method='sd') DataFrame

计算 eps 参数的收缩量。

eps 参数收缩量公式

若为标准差形式,计算公式如下:

\[ shr = [1 - SD(IWRES)] \times 100\% \]

若为方差形式,计算公式如下:

\[ shr = [1 - Var(IWRES)] \times 100\% \]

其中,IWRES 为个体加权残差(individual weighted residual)。

参数:

  • method (Literal["sd", "var"], optional):计算收缩量的方法。使用标准差形式("sd")或方差形式("var")。默认值为 "sd"。计算公式可见:eps 参数收缩量公式

返回值:

eps 参数收缩量的数据框。

返回类型:

DataFrame

5.3.19. plot_gof()#

plot_gof(self, filter=None, idv_name=None, dv_name=None, wres='CWRES', iwres='CIWRES') Plot#

绘制模型拟合优度图(goodness of fit)。

参数:

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

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

  • dv_name (str, optional):自变量名。若为 None,则自动推断。默认值为 None

  • wres (Literal["CWRES", "CWRESI", "WRES", "WRESI"], optional):需要计算的加权残差类型(weighted residuals)。默认值为 "CWRES"

  • iwres (Literal["CIWRES", "CIWRESI", "IWRES", "IWRESI"], optional):需要计算的个体加权残差类型(individual weighted residual)。默认值为 "CIWRES"

示例:

fit_res.plot_gof()

5.3.20. plot_predictions()#

plot_predictions(self, filter=None, idv_name=None, dv_name=None, wres='CWRES', iwres='CIWRES') Plot#

绘制个体拟合图(个体预测值/群体预测值 vs. 实测值)。

参数:

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

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

  • dv_name (str, optional):自变量名。若为 None,则自动推断。默认值为 None

  • wres (Literal["CWRES", "CWRESI", "WRES", "WRESI"], optional):需要计算的加权残差类型(weighted residuals)。默认值为 "CWRES"

  • iwres (Literal["CIWRES", "CIWRESI", "IWRES", "IWRESI"], optional):需要计算的个体加权残差类型(individual weighted residual)。默认值为 "CIWRES"

示例:

fit_res.plot_predictions(filter=pl.col("SEX") == 0)

5.3.21. plot_trace()#

plot_trace(self) Plot#

绘制迭代历史图(目标函数值与各模型参数估计值在迭代间的变动情况)。

示例:

fit_res.plot_trace()
c:\Users\user\workspace\mas-dev\mas\libs\phanpy\plotting1\scales\scales.py:60: PlotWarning: Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

5.3.22. vpc()#

vpc(n_sample=1000, seed=1234, ci_level=0.95, percentile_level=0.95, auto_bin=True, min_points_in_each_bin=10, linear_interpolation=True, *, dv_name=None, idv_name=None) VpcResult#

运行模型可视化预测检验(Visual Predictive Check)。

参数:

  • n_sample (int, optional):抽样次数。默认值为 1000

  • seed (int, optional):随机种子号。默认值为 1234

  • ci_level (float, optional):模拟数据百分位数的置信水平。默认值为 0.95

  • percentile_level (float, optional):模拟数据与实测值的百分位数水平。默认值为 0.95,即取 2.5% 与 97.5% 分位数。

  • auto_bin (bool | int | tuple[int, int], optional):分箱(binning)方法。默认值为 True。具体方法如下:

    • 若为 True,则使用自动分箱方法。

    • 若为 False,则视每一个观测时间点为一个分箱。

    • 若为 int,则视分箱数目为 int,将基于此数目计算每个箱体的范围。

    • 若为 tuple[int, int],则视箱体的范围为 tuple[int, int],将基于此范围计算最终的分箱数目。

  • min_points_in_each_bin (int, optional):每个分箱中需包含的最少观测点数目。默认值为 10

  • linear_interpolation (bool, optional):是否进行线性化。默认值为 True

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

  • dv_name (str, optional):自变量名。若为 None,则自动推断。默认值为 None

返回值:

包含 VPC 图、置信区间数据、模拟数据等相关内容的检验结果。

返回类型:

VpcResult

示例:

vpc_res = fit_res.vpc()
vpc_res.plot(linear_interpolation=False)
[VPC] Preparing data...
[VPC] Running simulation...
[VPC] Done in 2.980676s
Running Auto-bin...
Bin search ended, result is 10
Auto-bin completed, finally use 10 bins.

5.3.23. bootstrap()#

bootstrap(n_sample=1000, seed=1234, level=0.95) BootstrapResult#

运行自举法(bootstrap)。

参数:

  • n_sample (int, optional):抽样次数。默认值为 1000

  • seed (int, optional):随机种子号。默认值为 1234

  • level (float, optional):参数估计值的置信水平。默认值为 0.95

返回值:

包含估算成功率与抽样数据等相关内容的检验结果。

返回类型:

BootstrapResult

5.3.24. scm()#

scm(spec) ScmResult#

利用 SCM 法(stepwise covariate modeling)进行协变量筛选。

参数:

  • spec (MarkedScmSpecDefFunctionType | ScmSpec):SCM 配置。可详见 ScmSpec

返回值:

SCM 运行结果汇总。

返回类型:

ScmResult