2016-06-07 16 views
0

Ich versuche, ein Emax-Modell basierend auf den Daten und das Modell in diesem Video .. (ca. 40 Minuten in)Erstellen von Emax Modell in pymc3

https://www.youtube.com/watch?v=U9Nf-ZYHRQA&feature=youtu.be&list=PLvLDbH2lpyXNGV8mpBdF7EFK9LQJzGL-Y

hier mit pymc3 zu bauen, ist ein Screenshot zeigt das Modell ... enter image description here

Mein Code ist hier ...

pkpd_model = Model() 

with pkpd_model: 

# Hyperparameter Priors 
mu_e0 = Normal('mu_e0', mu=0, tau =1000) 
tau_e0 = Uniform('tau_e0', lower=0, upper =100) 

mu_emax = Normal('mu_emax', mu=0, tau =1000) 
tau_emax = Uniform('tau_emax', lower=0, upper =100) 

e0 = Lognormal('e0', mu = mu_e0, tau=tau_e0, shape =n_studies) 
emax= Lognormal('emax', mu = mu_emax, tau =tau_emax, shape =n_studies) 
ed50 = Lognormal('ed50', mu=1, tau = 1000) 

# Normalise sigma for sample size 
sigma = np.sqrt(np.square(Uniform('sigma', lower = 0, upper = 1000))/n) 


# Expected value of outcome 
resp_median = e0[study] + (emax[study]*dose)/(ed50+dose) 


# Likelihood (sampling distribution) of observations 
resp = Lognormal('resp', mu=resp_median, tau =sigma, observed =mean_response) 
resp_pred = Lognormal('resp_pred', mu=resp_median, tau =sigma, shape =len(dose)) 

das Modell r uns passt und passt, es ist nur so, dass die hinteren Schätzungen für den Modellparameter nicht in der Nähe sind, was ich erwarte. zum Beispiel ist meine Schätzung von emax ungefähr 2, aber Sie können klar sehen, dass von den Daten es ungefähr 10 sein sollte. So kann ich nur annehmen, dass ich einen Fehler gemacht habe, das Modell zu bauen, aber ich kann nicht für das Leben von mir sehen, was es ist.

Können Sie helfen?

Dank

Mark

Antwort

0

Nachdem die Modelldaten Neuüberprüfen sah ich, dass die resp_median tatsächlich das Protokoll dieses Ausdrucks ist. So ...

esp_median = np.log(e0[study] + (emax[study]*dose)/(ed50+dose)) 

verbessert Dinge. Ich bin mir aber immer noch nicht sicher, ob ich die Distributionen korrekt angegeben habe. Irgendwelche Hinweise irgendjemand?