Creating Emax model in pymc3

Question

I am trying to build an Emax model using pymc3 based on the data and model in this video.. (about 40mins in)

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

Here is a screen shot showing the model...

My code is here...

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

The model runs and fits okay, its just that the posterior estimates for the model parameter are not near what I am expecting. for example my emax estimate is around 2 but you can clearly see that from the data it should be around 10. So I can only assume I have made an error in building the model, but I cannot for the the life of me see what it is.

Can you help?

Thanks

Mark

Answer 1

After reexamining the model specs I saw that the resp_median is actually the log of that expression. So...

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

improves things. However I am still not sure that I have specified the distributions correctly. Any pointers anyone?