Ich bin neu in der PK-Modellierung und in pymc3, aber ich habe mit pymc3 herumgespielt und versucht, ein einfaches PK-Modell als Teil meines eigenen Lernens zu implementieren. Speziell ein Modell, das diese Beziehung einfängt ...PyMC3 PK Modellierung. Modell kann nicht zu Parametern aufgelöst werden, die zum Erstellen des Datensatzes verwendet wurden
Wo C (t) (Cpred) Konzentration zur Zeit t ist, Dosis ist die Dosis, gegeben, V Verteilungsvolumen ist, ist CL Clearance.
Ich habe einige Testdaten (30 Probanden) mit Werten von CL = 2, V = 10, für 3 Dosen 100,200,300 und generierte Daten zu den Zeitpunkten 0,1,2,4,8,12 generiert und auch enthalten ein zufälliger Fehler auf CL (Normalverteilung, 0 Mittel, Omega = 0,6) und auf dem verbleibenden unerklärten Fehler DV = Cpred + Sigma, wobei Sigma normalverteilt ist, SD = 0,33. Zusätzlich habe ich eine Transformation von C und V in Bezug auf das Gewicht (gleichmäßige Verteilung 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']]
Arrays für das Modell fertig erstellen ...
# 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
ich ein einfaches Modell in pymc3 gebaut haben ....
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)
Meine Erwartung war, dass ich in der Lage sein sollte um aus den Daten die Konstanten und Fehlerraten zu ermitteln, die ursprünglich zur Erzeugung der Daten verwendet wurden, obwohl ich dies nicht konnte, obwohl ich mich nähern kann. In diesem Beispiel ...
data['CLi'].mean()
> 2.322473543135788
data['Vi'].mean()
> 10.147619047619049
Und die Spur zeigt ....
Also meine Fragen sind ..
- Ist mein Code richtig strukturiert und gibt es irgendwelche eklatanten Fehler, die ich übersehen habe, könnten für diesen Unterschied verantwortlich sein?
- Kann ich das Modell pymc3 strukturieren, um die Beziehung, aus der ich die Daten generiert habe, besser widerzuspiegeln?
- Was wären Ihre Vorschläge zur Verbesserung des Modells?
Vielen Dank im Voraus!