PyMC3 PKモデリング。モデルはデータセットを作成するために使用されるパラメータに解決されません

Question

私はPKモデリングとpymc3の新人ですが、私はpymc3で遊んでいて、自分の学習の一部として単純なPKモデルを実装しようとしています。具体的には、この関係を捕捉するモデル... C（T）（Cpred）は時間tでの濃度である

、用量は、所与の用量であり、Vは分布容積であり、CLはクリアランスです。

私は、CL = 2、V = 10の値を3回の投与量100,200,300で生成し、時点0,1,2,4,8,12でデータを生成し、 （正規分布、0平均、オメガ= 0.6）および残存する原因不明の誤差DV = Cpred + sigma（ここで、シグマは通常SD = 0.33に分布している）に関するいくつかのランダム誤差。さらに、私は、重量（一様分布50-90）CLi = CL * WT/70に関してCとVの変換を含めた。 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']] 

私は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) 

私の期待は、私ができているはずということでした...

# 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 

をモデルの準備の配列を作成します。私はこれを行うことができませんでしたが、データを生成するためにもともと使用されていた定数とエラーレートをデータから解決しました。この例では...

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

とトレースショー....

だから私の質問は...

1. ある私のコードが正しく構成され、そこにあるされています私が見落とした目障りなミスは、この違いを説明してくれるでしょうか？
2. データを生成した関係をよりよく反映するようにpymc3モデルを構成できますか？
3. モデルを改善するための提案はありますか？

ありがとうございます！

Answer 1

私自身の質問にお答えします！

しかし、私はここで見つけるたとえば、次の階層モデルを実装し...

GLM -hierarchical

、それが御馳走に動作します。また、私は、私がデータフレーム内のエラーを適用した方法でエラーに気づいた - 各被験者を確保するために

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

を使用する必要がありますエラー

に異なる値を持つ