PyMC3 PK 建模。模型无法解析用于创建数据集的参数

Question

我是 PK 建模和 pymc3 的新手，但我一直在玩 pymc3 并尝试实现一个简单的 PK 模型作为我自己学习的一部分。特别是捕捉这种关系的模型...

其中 C(t)(Cpred) 是时间 t 的浓度，Dose 是给定的剂量，V 是分布容积，CL 是清除率。

我生成了一些测试数据（30 名受试者），CL =2，V=10，3 剂 100,200,300，并在时间点 0、1、2、4、8、12 生成了数据，并且在 CL（正态分布，0 均值，omega =0.6）和剩余无法解释的误差 DV = Cpred + sigma 上包括一些随机误差，其中 sigma 是正态分布的，SD =0.33。此外，我还对 C 和 V 进行了权重转换（均匀分布 50-90） CLi = CL * WT/70； Vi = V * WT/70.

# Create Data for modelling
np.random.seed(0)
# Subject ID's
data = pd.DataFrame(np.arange(1,31), columns=['subject'])
# Dose
Data['dose'] = np.array([100,100,100,100,100,100,100,100,100,100,
                    200,200,200,200,200,200,200,200,200,200,
                    300,300,300,300,300,300,300,300,300,300])
# Random Body Weight
data['WT'] = np.random.randint(50,100, size =30)
# Fixed Clearance and Volume for the population
data['CLpop'] =2
data['Vpop']=10
# Error rate for individual clearance rate
OMEGA = 0.66
# Individual clearance rate as a function of weight and omega
data['CLi'] = data['CLpop']*(data['WT']/70)+ np.random.normal(0, OMEGA )
# Individual Volume as a function of weight
data['Vi'] = data['Vpop']*(data['WT']/70) 

# Expand dataframe to account for time points
data = pd.concat([data]*6,ignore_index=True )
data = data.sort('subject')
# Add in time points
data['time'] = np.tile(np.array([0,1,2,4,8,12]), 30)

# Create concentration values using equation
data['Cpred'] = data['dose']/data['Vi'] *np.exp(-1*data['CLi']/data['Vi']*data['time'])
# Error rate for DV
SIGMA = 0.33
# Create Dependenet Variable from Cpred + error
data['DV']= data['Cpred'] + np.random.normal(0, SIGMA )

# Create new df with only data for modelling...
df = data[['subject','dose','WT', 'time', 'DV']]

为模型创建阵列...

# Prepare data from df to model specific arrays
time = np.array(df['time'])
dose = np.array(df['dose'])
DV = np.array(df['DV'])
WT = np.array(df['WT'])
n_patients = len(data['subject'].unique())
subject = data['subject'].values-1

我已经在 pymc3 中构建了一个简单的模型....

pk_model = Model()

with pk_model:
    # Hyperparameter Priors     
    sigma = Lognormal('sigma', mu =0, tau=0.01)

    V = Lognormal('V', mu =2, tau=0.01)
    CL = Lognormal('CL', mu =1, tau=0.01)    

    # Transformation wrt to weight
    CLi = CL*(WT)/70    
    Vi = V*(WT)/70

    # Expected value of outcome
    pred = dose/Vi*np.exp(-1*(CLi/Vi)*time)

    # Likelihood (sampling distribution) of observations
    conc = Normal('conc', mu =pred, tau=sigma, observed = DV)

我的期望是我应该能够从数据中解析出最初用于生成数据的常数和错误率，虽然我还没有做到这一点，尽管我可以接近。在这个例子中...

data['CLi'].mean()
> 2.322473543135788
data['Vi'].mean()
> 10.147619047619049

跟踪显示....

所以我的问题是..

我的代码结构是否正确，是否有任何我忽略的明显错误可能导致这种差异？
我能否构建 pymc3 模型以更好地反映我生成数据的关系？
您对改进模型有何建议？

提前致谢！

Answer 1

我要回答我自己的问题了！

但我按照此处找到的示例实现了分层模型...

GLM -hierarchical

而且效果很好。我还注意到我在数据框中应用错误的方式存在错误 - 应该使用

data['CLer'] = np.random.normal(scale=OMEGA, size=30)

确保每个主题都有不同的错误值

PyMC3 PK 建模。模型无法解析用于创建数据集的参数

PyMC3 PK modelling. Model cant resolve to parameters used to create the data set

bayesian

pymc

pymc3