#!/usr/bin/env python
# coding: utf-8

# # Infantile Hemangioma Meta-analysis
# 
# 
# > KQ 2: Among newborns, infants, and children up to 18 years of age with infantile hemangiomas who have been referred for pharmacologic intervention, what is the comparative effectiveness (benefits/harms) of corticosteroids or beta-blockers?
# 
# Using data extracted by the systematic review, we will conduct a meta-analysis to estimate the effectiveness of several corticosteroids and beta blockers for the treatment of infantile hemangiomas. Of particular interest is the estimation of the efficacy of propanolol, a beta blocker that was used in a large number of studies in the review. To this end, we will determine the expected probability of partial-to-complete clearance for propanolol, and compare this to the same probabilities for several different comparators, as determined by their availability in the dataset.
# 
# ### Beta-binomial model
# 
# For this model, the response variable is the number of individuals in study $j$ under intervention $k$ that achieve the clearance threshold:
# 
# $$y_{jk} = \sum_{i=1}^{n_{jk}} = I_i(\text{above clearance threshold})$$
# 
# This outcome is modeled as a binomial response:
# 
# $$y_{jk} \sim \text{Bin}(n_{jk}, \pi_{jk})$$
# 
# where $\pi_{jk})$ is the probability of a positive response for study $j$ under intervention $k$. To allow for heterogeneity in this probability across studies, we can specify it as a random effect:
# 
# $$\pi_{jk} \sim \text{Beta}(\alpha, \beta)$$
# 
# where $\alpha, \beta$ are the parameters of a beta distribution (which models quantities on the [0,1] interval), resulting in a beta-binomial distribution for the outcome.
# 
# It may be possible to incorporate covariates to improve the prediction of propanolol effectiveness. In particular, the mode of delivery (oral, intralesional, topical), dose, or the hemangioma location may be predictive of intervention effectiveness. For a vector of such covariates $x$, we can alternatively model $\pi_{jk}$ as a logit-linear function:
# 
# $$\text{logit}(\pi_{jk}) = x_j^{\prime}\theta_k + \epsilon_j$$
# 
# where $\theta$ is a corresponding vector of regression parameters corresponding to intervention $k$, and $\epsilon_j$ is a study-level random effect to account for the correlation of study arms.
# 
# ### Latent variable model
# 
# However, the use of an arbitrary cutoff value as a threshold of success is an unfortunate and perhaps unsatisfactory modeling choice. There is an inherent loss of information in the dichotomization of continuous variables, and this loss is magnified here by having to discard data from studies that use a different response threshold than the adopted value (*e.g.* 75%). Since the clearance rate is a continuous measure, one can hypothesize a latent, continuous probability distribution that each study reports according to its respective quantiles: 25%, 50%, 75%, etc. If there is sufficient information, one may use a Bayesian approach to attempt to reconstruct this latent distribution, which would allow for more of the available information to be used in the meta-analytic procedure.  
# 
# Under some treatment $k$, one can consider a notional *distribution* of hemangioma clearance rates, from no effect (0) to complete clearance (1). As a matter of convenience in any particular study $j$, researchers will have chosen an arbitrary clearance threshold $c_j$, only recording whether a particular subject occupied one side or the other of this threshold. We can characterize the true, latent response distribution by estimating the parameters $\mu_k, \sigma_k$ via the following identity:
# 
# $$\pi_{jk} = 1 - I_{c_j}(\mu_k, \sigma_k)$$
# 
# where $I_x(a,b)$ is the cumulative distribution function of a logit-transformed normal distribution under parameters $\mu, \sigma$. The resulting probability is the same as specified above, and can be used in the same binomial likelihood:
# 
# $$y_{jk} \sim \text{Bin}(n_{jk}, \pi_{jk})$$
# 
# ### Comparative effectiveness
# 
# Irrespective of which model form is employed, the comparative effectiveness of any two pharmacologic interventions $k=1,2$ can be assessed directly by the difference in their respective clearance threshold probabilities:
# 
# $$d_{12} = \pi_1 - \pi_2$$
# 
# Using a Bayesian framework, we can extract associated 95% posterior credible intervals for this difference.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sb
import pylab as pl

from scipy.stats import norm


# In[132]:


sb.set_style("white")


# In[2]:


# Set random number seed
np.random.seed(42)


# In[3]:


kq2_outcomes = pd.read_excel('IH_effectiveness_data_Extraction_comparative studies_revised_03242015.xlsx', 
                             sheetname='KQ2_outcome data')

# Individual patient data for RefID 112
kq2_ipd = pd.read_excel('IH_effectiveness_data_Extraction_comparative studies_revised_03242015.xlsx', 
                             sheetname='KQ2_IPD')


# Individual patient data for REFID 112 (Qiu 2013). Individual with multiple treatments per lesion was deleted; data from multiple lesions on the same individual were averaged.

# In[4]:


unique_outcomes = kq2_outcomes['Outcome'].unique()


# In[5]:


unique_interventions = kq2_outcomes['Arm'].unique()

unique_interventions


# In[6]:


intervention_categories = pd.read_csv('interventions.csv', index_col=0)


# Remove unusable interventions

# In[7]:


dropped_interventions = intervention_categories[intervention_categories.Label=='DROP'].index.values
kq2_outcomes = kq2_outcomes[~kq2_outcomes.Arm.isin(dropped_interventions)]


# Recode interventions into collapsed category set

# In[8]:


intervention_lookup = intervention_categories[intervention_categories.Label!='DROP'].to_dict()['Label']


# In[9]:


intervention_lookup


# In[10]:


kq2_outcomes['intervention'] = kq2_outcomes.Arm.replace(intervention_lookup)
kq2_outcomes['intervention'].value_counts()


# In[11]:


unique_interventions = kq2_outcomes.intervention.unique()
n_interventions = len(unique_interventions) - 1
n_interventions


# In[12]:


unique_interventions


# Define indices to intervention elements

# In[13]:


PROPRANOLOL = 0
ATENOLOL = 1
TIMOLOL = 2
TRIAMCILONE = 3
ORAL_STEROID = 4
IMIQUIMOD = 5
NADOLOL = 6
METHYLPREDNISOLONE = 7


# Number of unique studies

# In[14]:


unique_studies = kq2_outcomes.REFID.unique()
n_studies = len(unique_studies)
n_studies


# In[15]:


unique_studies


# There does not appear to be sufficient variation in dosage to estimate its effect for any agent.

# In[16]:


kq2_outcomes.groupby(['Dosage form','Dosage amount','intervention'])['Dosage'].value_counts()


# In[17]:


kq2_outcomes.groupby(['intervention', 'Dosage form']).REFID.count()


# ## Model specification
# 
# Study random effect. Might want to use a robust distribution, like Cauchy or t.

# In[18]:


σ_ϵ = pm.Uniform('σ_ϵ', 0, 1000, value=10)
τ_ϵ = σ_ϵ ** -2

ϵ = pm.Normal('ϵ', 0, τ_ϵ, value=np.zeros(n_studies))


# Parameters for each intervention

# In[19]:


μ = pm.Normal('μ', 0, 0.001, value=-2)
θ = pm.Normal('θ', 0, 0.001, value=np.zeros(n_interventions))
λ = pm.Exponential('λ', 0.1, value=0.5)
σ = pm.HalfNormal('σ', λ, value=np.ones(n_interventions+1))
τ = σ ** -2


# In[20]:


θ.value.shape


# In[21]:


σ.value.shape


# Alternate mode of delivery effects for propanolol (oral is baseline)

# In[22]:


ϕ = pm.Normal('ϕ', 0, 0.001, value=0)


# ### Abarrzua-Araya 2014 (RefID 5)

# In[23]:


subset_cols = ['REFID', 'Arm', 'Dosage', 'Dosage amount', 'Dosage form', 'N at follow-up', 'Outcome', 
    'Outcome data/result', 'Outcome data', 'Point estimate', 'sd', 'Lower range', 'Upper range', 'Outcome data/result N']


# In[24]:


kq2_outcomes[kq2_outcomes.REFID==5][subset_cols]


# In[25]:


refid5 = np.where(unique_studies==5)[0][0]


# In[26]:


@pm.deterministic
def mu_propranolol_5(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for propanolol in RefID 5
    return μ + θ[PROPRANOLOL] + ϵ[refid5]


@pm.deterministic
def p_propranolol_5(μ=mu_propranolol_5, τ=τ):
    # Response category probabilities for propanolol in RefID 5
    c1 = norm.cdf(pm.logit(0.01), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.99), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, 1-c2   

# Multinomial likelihood for propanolol in RefID 5
y_propanolol_5 = pm.Multinomial('y_propanolol_5', 10, p_propranolol_5, value=[0, 4, 6], observed=True)


# In[27]:


@pm.deterministic
def mu_atenolol_5(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for atenolol in RefID 5
    return μ + θ[ATENOLOL] + ϵ[refid5]


@pm.deterministic
def p_atenolol_5(μ=mu_atenolol_5, τ=τ):    
    # Response category probabilities for atenolol in RefID 5
    c1 = norm.cdf(pm.logit(0.01), μ, τ[ATENOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.99), μ, τ[ATENOLOL]**-1)
    return c1, c2-c1, 1-c2  

# Multinomial likelihood for atenolol in RefID 5
y_atenolol_5 = pm.Multinomial('y_atenolol_5', 13, p_atenolol_5, value=[0, 6, 7], observed=True)


# ### Chan 2013 (RefID 110)

# In[28]:


kq2_outcomes[kq2_outcomes.REFID==110][subset_cols]


# In[29]:


refid110 = np.where(unique_studies==110)[0][0]


# In[30]:


@pm.deterministic
def mu_timolol_110(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for timolol in RefID 110
    return μ + θ[TIMOLOL] + ϵ[refid110]


@pm.deterministic
def p_timolol_110(μ=mu_timolol_110, τ=τ):
    # Response category probabilities for timolol in RefID 110
    return 1 - norm.cdf(pm.logit(0.05), μ, τ[TIMOLOL]**-1)

# Binomial likelihood for propanolol in RefID 5
y_timolol_110 = pm.Binomial('y_timolol_110', 15, p_timolol_110, value=9, observed=True)


# In[31]:


@pm.deterministic
def mu_control_110(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for control in RefID 110
    return μ + ϵ[refid110]


@pm.deterministic
def p_control_110(μ=mu_control_110, τ=τ):
    # Response category probabilities for control in RefID 110
    return 1 - norm.cdf(pm.logit(0.05), μ, τ[-1]**-1)

# Binomial likelihood for control in RefID 5
y_control_110 = pm.Binomial('y_control_110', 18, p_control_110, value=2, observed=True)


# ### Jalil 2006 (RefID 1263)

# In[32]:


kq2_outcomes[kq2_outcomes.REFID==1263][subset_cols]


# In[33]:


refid1263 = np.where(unique_studies==1263)[0][0]


# In[34]:


@pm.deterministic
def mu_triamcilone_1263(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for triamcilone in RefID 1263
    return μ + θ[TRIAMCILONE] + ϵ[refid1263]


@pm.deterministic
def p_triamcilone_1263(μ=mu_triamcilone_1263, τ=τ):
    # Response category probabilities for triamcilone in RefID 1263
    c1 = norm.cdf(pm.logit(0.1), μ, τ[TRIAMCILONE]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[TRIAMCILONE]**-1)
    return c1, c2-c1, 1-c2   

# Multinomial likelihood for triamcilone in RefID 1263
y_triamcilone_1263 = pm.Multinomial('y_triamcilone_1263', 25, p_triamcilone_1263, value=[6, 8, 11], observed=True)


# In[35]:


@pm.deterministic
def mu_prednisolone_1263(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for prednisolone in RefID 1263
    return μ + θ[ORAL_STEROID] + ϵ[refid1263]


@pm.deterministic
def p_prednisolone_1263(μ=mu_triamcilone_1263, τ=τ):
    # Response category probabilities for prednisolone in RefID 1263
    c1 = norm.cdf(pm.logit(0.1), μ, τ[ORAL_STEROID]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[ORAL_STEROID]**-1)
    return c1, c2-c1, 1-c2   

# Multinomial likelihood for prednisolone in RefID 1263
y_prednisolone_1263 = pm.Multinomial('y_prednisolone_1263', 25, p_prednisolone_1263, value=[6, 11, 8], observed=True)


# In[36]:


@pm.deterministic
def mu_control_1263(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for control in RefID 1263
    return μ + ϵ[refid1263]


@pm.deterministic
def p_control_1263(μ=mu_control_1263, τ=τ):
    # Response category probabilities for prednisolone in RefID 1263
    c1 = norm.cdf(pm.logit(0.01), μ, τ[-1]**-1)
    c2 = norm.cdf(pm.logit(0.1), μ, τ[-1]**-1)
    c3 = norm.cdf(pm.logit(0.5), μ, τ[-1]**-1)
    return c1, c2-c1, c3-c2, 1-c3 

# Multinomial likelihood for control in RefID 1263
y_control_1263 = pm.Multinomial('y_control_1263', 25, p_control_1263, value=[5, 19, 1, 0], observed=True)


# ### Zahar 2013 (RefID 43)

# In[37]:


kq2_outcomes[kq2_outcomes.REFID==43][subset_cols]


# In[38]:


refid43 = np.where(unique_studies==43)[0][0]


# In[39]:


@pm.deterministic
def mu_propranolol_oral_43(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 43
    return μ + θ[PROPRANOLOL] + ϵ[refid43]

@pm.deterministic
def mu_propranolol_oral_43(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 43
    return μ + θ[PROPRANOLOL] + ϵ[refid43]
@pm.deterministic
def p_propranolol_oral_43(μ=mu_propranolol_oral_43, τ=τ):
    # Response category probabilities for oral propanolol in RefID 43
    c1 = norm.cdf(pm.logit(0.01), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.50), μ, τ[PROPRANOLOL]**-1)
    c3 = norm.cdf(pm.logit(0.99), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, c3-c2, 1-c3  

# Multinomial likelihood for oral propanolol in RefID 43
y_propanolol_oral_43 = pm.Multinomial('y_propanolol_oral_43', 15, p_propranolol_oral_43, value=[3, 1, 2, 9], 
                                      observed=True)


# In[40]:


# @pm.deterministic
# def mu_propranolol_topical_43(μ=μ, θ=θ, ϕ=ϕ, ϵ=ϵ):
#     # Mean response on the logit scale for topical propanolol in RefID 43
#     return μ + θ[PROPRANOLOL] + ϕ + ϵ[refid43]


# @pm.deterministic
# def p_propranolol_topical_43(μ=mu_propranolol_topical_43, τ=τ):
#     # Response category probabilities for topical propanolol in RefID 43
#     c1 = norm.cdf(pm.logit(0.01), μ, τ**-1) #[PROPRANOLOL]**-1)
#     c2 = norm.cdf(pm.logit(0.50), μ, τ**-1) #[PROPRANOLOL]**-1)
#     c3 = norm.cdf(pm.logit(0.99), μ, τ**-1) #[PROPRANOLOL]**-1)
#     return c1, c2-c1, c3-c2, 1-c3  

# # Multinomial likelihood for topical propanolol in RefID 43
# y_propanolol_topical_43 = pm.Multinomial('y_propanolol_topical_43', 15, p_propranolol_topical_43, value=[5, 2, 5, 3], 
#                                       observed=True)


# In[41]:


@pm.deterministic
def mu_propranolol_intralesional_43(μ=μ, θ=θ, ϕ=ϕ, ϵ=ϵ):
    # Mean response on the logit scale for intralesional propanolol in RefID 43
    return μ + θ[PROPRANOLOL] + ϕ + ϵ[refid43]


@pm.deterministic
def p_propranolol_intralesional_43(μ=mu_propranolol_intralesional_43, τ=τ):
    # Response category probabilities for intralesional propanolol in RefID 43
    c1 = norm.cdf(pm.logit(0.01), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.50), μ, τ[PROPRANOLOL]**-1)
    c3 = norm.cdf(pm.logit(0.99), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, c3-c2, 1-c3  

# Multinomial likelihood for intralesional propanolol in RefID 43
y_propanolol_intralesional_43 = pm.Multinomial('y_propanolol_intralesional_43', 15, p_propranolol_intralesional_43, 
                                               value=[7, 3, 3, 2], observed=True)


# ### Bauman 2014 (RefID 13)
# 
# This data was given as a mean change on TSA, with 95% confidence interval. We will use a normal factor potential on augmented data based on the estimated mean, and the back-calculated standard deviation.

# In[42]:


kq2_outcomes[kq2_outcomes.REFID==13][subset_cols]


# In[43]:


refid13 = np.where(unique_studies==13)[0][0]


# In[44]:


_propranolol_sd = 9 * (0.83 - 0.44) / 3.92
_prednisolone_sd = 6 * (0.72 - 0.10) / 3.92


# In[45]:


propranolol_13_var = (_propranolol_sd**2) / (0.64 * (1-0.64))**2
prednisolone_13_var = (_prednisolone_sd**2) / (0.41 * (1-0.41))**2


# In[46]:


@pm.deterministic
def mu_propranolol_13(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 13
    return μ + θ[PROPRANOLOL] + ϵ[refid13]

propranolol_obs_13 = pm.Normal('propranolol_obs_13', mu_propranolol_13, propranolol_13_var**-1, 
                               value=pm.logit(0.64), observed=True)


# In[47]:


@pm.deterministic
def mu_prednisolone_13(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for prednisolone in RefID 43
    return μ + θ[ORAL_STEROID] + ϵ[refid13]

prednisolone_obs_13 = pm.Normal('prednisolone_obs_13', mu_prednisolone_13, prednisolone_13_var**-1, 
                                value=pm.logit(0.41), observed=True)


# ### Hogeling 2011 (RefID 445)

# In[48]:


kq2_outcomes[kq2_outcomes.REFID==445][subset_cols]


# Hogeling et al did not report a measure of uncertainty with the effect sizes, therefore this study must be excluded from the meta-analysis.

# ### Leaute-Labreze 2015 (RefID 3740)

# In[49]:


kq2_outcomes[kq2_outcomes.REFID==3740][subset_cols]


# In[50]:


refid3740 = np.where(unique_studies==3740)[0][0]


# In[51]:


@pm.deterministic
def mu_propanolol_3740(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for propanolol in RefID 3740
    return μ + θ[PROPRANOLOL] + ϵ[refid3740]


@pm.deterministic
def p_propanolol_3740(μ=mu_propanolol_3740, τ=τ):
    # Response category probabilities for propanolol in RefID 3740
    return 1 - norm.cdf(pm.logit(0.9), μ, τ[PROPRANOLOL]**-1)

# Binomial likelihood for propanolol in RefID 5
y_propanolol_3740 = pm.Binomial('y_propanolol_3740', 101, p_propanolol_3740, value=61, 
                                observed=True)


# In[52]:


@pm.deterministic
def mu_control_3740(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for controls in RefID 3740
    return μ + ϵ[refid3740]


@pm.deterministic
def p_control_3740(μ=mu_control_3740, τ=τ):
    # Response category probabilities for controls in RefID 3740
    return 1 - norm.cdf(pm.logit(0.9), μ, τ[-1]**-1)

# Binomial likelihood for propanolol in RefID 5
y_control_3740 = pm.Binomial('y_control_3740', 25, p_control_3740, value=2, 
                                observed=True)


# ### Sondhi 2013 (RefID 72)
# 
# The paper claims `n=14` for controls, but only reports 12 outcomes.

# In[53]:


kq2_outcomes[kq2_outcomes.REFID==72][subset_cols]


# In[54]:


refid72 = np.where(unique_studies==72)[0][0]


# In[55]:


@pm.deterministic
def mu_propranolol_72(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for propanolol in RefID 72
    return μ + θ[PROPRANOLOL] + ϵ[refid72]


@pm.deterministic
def p_propranolol_72(μ=mu_propranolol_72, τ=τ):
    # Response category probabilities for propanolol in RefID 72
    c1 = norm.cdf(pm.logit(0.1), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, 1-c2   

# Multinomial likelihood for propanolol in RefID 5
y_propanolol_72 = pm.Multinomial('y_propanolol_72', 31, p_propranolol_72, value=[3, 0, 28], 
                                observed=True)


# In[56]:


@pm.deterministic
def mu_control_72(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for controls in RefID 72
    return μ + ϵ[refid72]


@pm.deterministic
def p_control_72(μ=mu_control_72, τ=τ):
    # Response category probabilities for controls in RefID 72
    c1 = norm.cdf(pm.logit(0.1), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, 1-c2   

# Multinomial likelihood for propanolol in RefID 5
y_control_72 = pm.Multinomial('y_control_72', 12, p_control_72, value=[6, 2, 4], 
                                observed=True)


# ### Qiu 2013 (RefID 112)
# 
# Qiu et al. report individual patient data, which we incorporate directly into the meta-analysis.

# In[57]:


kq2_outcomes[kq2_outcomes.REFID==112][subset_cols]


# In[58]:


refid112 = np.where(unique_studies==112)[0][0]


# In[59]:


obs_imiquimod_112 = kq2_ipd[kq2_ipd.Arm=='imiquimod'].Outcome.values/100.
obs_timolol_112 = kq2_ipd[kq2_ipd.Arm=='timolol'].Outcome.values/100.


# In[60]:


kq2_ipd[kq2_ipd.Arm=='imiquimod'].Outcome.median()


# In[61]:


@pm.deterministic
def mu_timolol_112(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for timolol in RefID 112
    return μ + θ[TIMOLOL] + ϵ[refid112]


y_timolol_112 = pm.Normal('y_timolol_112', mu_timolol_112, τ[TIMOLOL]**-1, 
                         value=pm.logit(obs_timolol_112 - 1e-6), observed=True)

@pm.deterministic
def mu_imiquimod_112(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for imiquimod in RefID 112
    return μ + θ[IMIQUIMOD] + ϵ[refid112]


y_imiquimod_112 = pm.Normal('y_imiquimod_112', mu_imiquimod_112, τ[IMIQUIMOD]**-1, 
                         value=pm.logit(obs_imiquimod_112), observed=True)


# ### Thayal 2012 (RefID 309)

# In[62]:


kq2_outcomes[kq2_outcomes.REFID==309][subset_cols]


# In[63]:


refid309 = np.where(unique_studies==309)[0][0]


# In[64]:


@pm.deterministic
def mu_propranolol_309(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for propanolol in RefID 309
    return μ + θ[PROPRANOLOL] + ϵ[refid309]


@pm.deterministic
def p_propranolol_309(μ=mu_propranolol_309, τ=τ):
    # Response category probabilities for propanolol in RefID 309
    c1 = norm.cdf(pm.logit(0.25), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[PROPRANOLOL]**-1)
    c3 = norm.cdf(pm.logit(0.75), μ, τ[PROPRANOLOL]**-1)
    c4 = norm.cdf(pm.logit(0.9), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, c3-c2, c4-c3, 1-c4  

# Multinomial likelihood for propanolol in RefID 309
y_propanolol_309 = pm.Multinomial('y_propanolol_309', 10, p_propranolol_309, 
                                  value=[1, 0, 3, 4, 2], observed=True)


# ### Chambers 2012 (RefID 321)

# In[65]:


kq2_outcomes[kq2_outcomes.REFID==321][subset_cols]


# In[66]:


refid321 = np.where(unique_studies==321)[0][0]


# In[67]:


@pm.deterministic
def mu_timolol_321(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for timolol in RefID 321
    return μ + θ[TIMOLOL] + ϵ[refid321]


@pm.deterministic
def p_timolol_321(μ=mu_timolol_321, τ=τ):
    # Response category probabilities for timolol in RefID 321
    c1 = norm.cdf(pm.logit(0.01), μ, τ[TIMOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[TIMOLOL]**-1)
    return c1, c2-c1, 1 - c2

# Binomial likelihood for propanolol in RefID 321
y_timolol_321 = pm.Multinomial('y_timolol_321', 13, p_timolol_321, value=[1, 4, 8], observed=True)


# In[68]:


@pm.deterministic
def mu_control_321(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for controls in RefID 321
    return μ + ϵ[refid321]


@pm.deterministic
def p_control_321(μ=mu_control_321, τ=τ):
    # Response category probabilities for controls in RefID 321
    c1 = norm.cdf(pm.logit(0.01), μ, τ[-1]**-1)
    c2 = norm.cdf(pm.logit(0.5), μ, τ[-1]**-1)
    return c1, c2-c1, 1 - c2

# Binomial likelihood for controls in RefID 321
y_control_321 = pm.Multinomial('y_control_321', 10, p_control_321, value=[9, 1, 0], observed=True)


# ### Bertrand 2011 (RefID 402)

# In[69]:


kq2_outcomes[kq2_outcomes.REFID==402][subset_cols]


# In[70]:


refid402 = np.where(unique_studies==402)[0][0]


# Find SD on the logit scale using the delta method

# In[71]:


propanolol_402_var = (0.2247**2) / (0.7873 * (1-0.7873))**2
prednisone_402_var = (0.1221**2) / (0.4482 * (1-0.4482))**2


# In[72]:


@pm.deterministic
def mu_propranolol_402(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 402
    return μ + θ[PROPRANOLOL] + ϵ[refid402]

propranolol_obs_402 = pm.Normal('propranolol_obs_402', mu_propranolol_402, propanolol_402_var**-1, 
                                value=pm.logit(0.7873), observed=True)


# In[73]:


@pm.deterministic
def mu_prednisolone_402(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for prednisolone in RefID 402
    return μ + θ[ORAL_STEROID] + ϵ[refid402]

prednisolone_obs_402 = pm.Normal('prednisolone_obs_402', mu_prednisolone_402, prednisone_402_var**-1, 
                                 pm.logit(0.4482), observed=True)


# ### Price 2011 (RefID 438)

# In[74]:


kq2_outcomes[kq2_outcomes.REFID==438][subset_cols]


# In[75]:


refid438 = np.where(unique_studies==438)[0][0]


# In[76]:


@pm.deterministic
def mu_propanolol_438(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for propanolol in RefID 438
    return μ + θ[PROPRANOLOL] + ϵ[refid438]


@pm.deterministic
def p_propranolol_438(μ=mu_propanolol_438, τ=τ):
    # Response category probabilities for propanolol in RefID 438
    return 1 - norm.cdf(pm.logit(0.75), μ, τ[PROPRANOLOL]**-1)

# Binomial likelihood for propanolol in RefID 5
y_propranolol_438 = pm.Binomial('y_propranolol_438', 68, p_propranolol_438, 
                                value=56, observed=True)


# In[77]:


@pm.deterministic
def mu_corticosteriod_438(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for corticosteriods in RefID 438
    return μ + θ[ORAL_STEROID] + ϵ[refid438]


@pm.deterministic
def p_corticosteriod_438(μ=mu_corticosteriod_438, τ=τ):
    # Response category probabilities for corticosteriods in RefID 438
    return 1 - norm.cdf(pm.logit(0.75), μ, τ[ORAL_STEROID]**-1)

# Binomial likelihood for propanolol in RefID 5
y_corticosteriod_438 = pm.Binomial('y_corticosteriod_438', 42, p_corticosteriod_438, 
                                value=12, observed=True)


# ### Yu 2013 (RefID 3451)

# In[78]:


kq2_outcomes[kq2_outcomes.REFID==3451][subset_cols]


# In[79]:


refid3451 = np.where(unique_studies==3451)[0][0]


# In[80]:


@pm.deterministic
def mu_timolol_3451(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for timolol in RefID 3451
    return μ + θ[TIMOLOL] + ϵ[refid3451]


@pm.deterministic
def p_timolol_3451(μ=mu_timolol_3451, τ=τ):
    # Response category probabilities for timolol in RefID 3451
    c1 = norm.cdf(pm.logit(0.01), μ, τ[TIMOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.99), μ, τ[TIMOLOL]**-1)
    return c1, c2-c1, 1 - c2

# Binomial likelihood for propanolol in RefID 3451
y_timolol_3451 = pm.Multinomial('y_timolol_3451', 101, p_timolol_3451, value=[44, 57-12, 12], observed=True)


# In[81]:


@pm.deterministic
def mu_control_3451(μ=μ, ϵ=ϵ):
    # Mean response on the logit scale for controls in RefID 3451
    return μ + ϵ[refid3451]


@pm.deterministic
def p_control_3451(μ=mu_control_3451, τ=τ):
    # Response category probabilities for timolol in RefID 3451
    c1 = norm.cdf(pm.logit(0.01), μ, τ[-1]**-1)
    c2 = norm.cdf(pm.logit(0.99), μ, τ[-1]**-1)
    return c1, c2-c1, 1 - c2

# Binomial likelihood for propanolol in RefID 3451
y_control_3451 = pm.Multinomial('y_control_3451', 23, p_control_3451, value=[22, 1, 0], observed=True)


# ### Awadein 2011 (RefID 3522)

# In[82]:


kq2_outcomes[kq2_outcomes.REFID==3522][subset_cols]


# In[83]:


refid3522 = np.where(unique_studies==3522)[0][0]


# In[84]:


@pm.deterministic
def mu_propranolol_intralesional_3522(μ=μ, θ=θ, ϕ=ϕ, ϵ=ϵ):
    # Mean response on the logit scale for intralesional propanolol in RefID 3522
    return μ + θ[PROPRANOLOL] + ϕ + ϵ[refid3522]


@pm.deterministic
def p_propranolol_intralesional_3522(μ=mu_propranolol_intralesional_3522, τ=τ):
    # Response category probabilities for intralesional propanolol in RefID 3522
    c1 = norm.cdf(pm.logit(0.01), μ, τ[PROPRANOLOL]**-1)
    c2 = norm.cdf(pm.logit(0.50), μ, τ[PROPRANOLOL]**-1)
    c3 = norm.cdf(pm.logit(0.99), μ, τ[PROPRANOLOL]**-1)
    return c1, c2-c1, c3-c2, 1-c3 

# Multinomial likelihood for intralesional propanolol in RefID 3522
y_propanolol_intralesional_3522 = pm.Multinomial('y_propanolol_intralesional_3522', 12, p_propranolol_intralesional_43, 
                                               value=[2, 2, 3, 5], observed=True)


# In[85]:


@pm.deterministic
def mu_triamcinolone_3522(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for triamcinolone in RefID 3522
    return μ + θ[TRIAMCILONE] + ϵ[refid3522]


@pm.deterministic
def p_triamcinolone_3522(μ=mu_triamcinolone_3522, τ=τ):
    # Response category probabilities for triamcinolone in RefID 3522
    c1 = norm.cdf(pm.logit(0.01), μ, τ[TRIAMCILONE]**-1)
    c2 = norm.cdf(pm.logit(0.50), μ, τ[TRIAMCILONE]**-1)
    c3 = norm.cdf(pm.logit(0.99), μ, τ[TRIAMCILONE]**-1)
    return c1, c2-c1, c3-c2, 1-c3 

# Multinomial likelihood for triamcinolone in RefID 3522
y_triamcinolone_3522 = pm.Multinomial('y_triamcinolone_3522', 10, p_triamcinolone_3522, 
                                               value=[2, 2, 2, 4], observed=True)


# ### Pope 2012 (RefID 3723)

# In[86]:


kq2_outcomes[kq2_outcomes.REFID==3723][subset_cols]


# In[87]:


refid3723 = np.where(unique_studies==3723)[0][0]


# In[88]:


propanolol_3723_var = (0.1482**2) / (0.86 * (1-0.86))**2
nadolol_3723_var = (0.0305**2) / (0.97 * (1-0.97))**2


# In[89]:


@pm.deterministic
def mu_propranolol_3723(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 402
    return μ + θ[PROPRANOLOL] + ϵ[refid3723]

propranolol_obs_3723 = pm.Normal('propranolol_obs_3723', mu_propranolol_3723, propanolol_3723_var**-1, 
                                value=pm.logit(0.86), observed=True)


# In[90]:


@pm.deterministic
def mu_nadolol_3723(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 402
    return μ + θ[NADOLOL] + ϵ[refid3723]

nadolol_obs_3723 = pm.Normal('nadolol_obs_3723', mu_nadolol_3723, nadolol_3723_var**-1, 
                                 value=pm.logit(0.97), observed=True)


# ### Pope 2007 (RefID 1180)

# In[91]:


kq2_outcomes[kq2_outcomes.REFID==1180][subset_cols]


# In[92]:


refid1180 = np.where(unique_studies==1180)[0][0]


# Assuming the width of the interquartile range is 1.35 SD.

# In[93]:


prednisolone_sd_1180 = (0.67 - 0.35) / 1.35
methylprednisolone_sd_1180 = (0.22 + 0.35) / 1.35


# In[94]:


methylprednisolone_1180_latent = pm.Normal('methylprednisolone_1180_latent', -0.015, 
                                             methylprednisolone_sd_1180**-2, size=10)


# In[95]:


methylprednisolone_1180_trunc = pm.Lambda('methylprednisolone_1180_trunc', 
                                         lambda y=methylprednisolone_1180_latent: np.clip(y, 1e-6, 1-1e-6))


# In[96]:


prednisolone_1180_latent = pm.Normal('prednisolone_1180_latent', 0.5, 
                                             prednisolone_sd_1180**-2, size=10)


# In[97]:


prednisolone_1180_trunc = pm.Lambda('prednisolone_1180_trunc', 
                                         lambda y=prednisolone_1180_latent: np.clip(y, 1e-6, 1-1e-6)).value


# In[98]:


@pm.deterministic
def mu_prednisolone_1180(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for prednisolone in RefID 1180
    return μ + θ[ORAL_STEROID] + ϵ[refid1180]

@pm.potential
def prednisolone_obs_1180(μ=mu_prednisolone_1180, τ=τ, y=prednisolone_1180_trunc):
    return pm.normal_like(pm.logit(y), μ, τ[ORAL_STEROID])


# In[99]:


@pm.deterministic
def mu_methylprednisolone_1180(μ=μ, θ=θ, ϵ=ϵ):
    # Mean response on the logit scale for oral propanolol in RefID 1180
    return μ + θ[METHYLPREDNISOLONE] + ϵ[refid1180]

@pm.potential
def methylprednisolone_obs_1180(μ=mu_methylprednisolone_1180, τ=τ, y=methylprednisolone_1180_trunc):
    return pm.normal_like(pm.logit(y), μ, τ[METHYLPREDNISOLONE])


# Expected clearance rates for each treatment

# In[100]:


@pm.deterministic
def expected_clearance(μ=μ, θ=θ, ϕ=ϕ):
    p = np.append(pm.invlogit(μ), pm.invlogit(μ + θ[0] + ϕ))
    return np.append(p, pm.invlogit(μ + θ))


# In[101]:


@pm.deterministic
def predicted_clearance(μ=μ, θ=θ, ϕ=ϕ, τ=τ):
    m = np.append(np.append(μ, μ + θ[0] + ϕ), μ + θ)
    T = np.append(np.append(τ[-1], τ[0]), τ[:-1])    
    y = pm.rnormal(m, T)
    return pm.invlogit(y)


# In[102]:


best = pm.Lambda('best', lambda b=expected_clearance:  (b==b.max()).astype(int))


# In[103]:


rate_labels = ['Control', 
               'Intralesional propanolol', 
               'Oral propanolol', 
               'Atenolol', 
               'Timolol', 
               'Triamcilone', 
               'Oral steroid',
               'Imiquimod', 
               'Nadolol', 
               'Methylprednisolone']


# In[ ]:


iterations, burn = 200000, 190000

M = pm.MCMC(locals())
M.use_step_method(pm.AdaptiveMetropolis, σ)
M.sample(iterations, burn)


# In[113]:


M.sample(iterations, burn)


# In[114]:


best.summary()


# In[125]:


pm.Matplot.plot(θ)


# In[115]:


pm.Matplot.plot(σ)


# In[116]:


μ.summary()


# In[117]:


θ.summary()


# In[118]:


pl.hist(pm.invlogit(μ.trace() + θ.trace()[:, 5]))


# In[136]:


sigma_labels = ['Propanolol', 
               'Atenolol', 
               'Timolol', 
               'Triamcilone', 
               'Oral steroid',
               'Imiquimod', 
               'Nadolol', 
               'Methylprednisolone',
               'Control']

pm.Matplot.summary_plot(σ, custom_labels=sigma_labels, 
                        vline_pos=-1, main='Standard deviation of effects')


# In[138]:


pm.Matplot.summary_plot(θ, custom_labels=sigma_labels, 
                        vline_pos=-10, main='Expected effect size (relative to control)')


# In[121]:


pm.Matplot.summary_plot(predicted_clearance, custom_labels=rate_labels, vline_pos=-1)


# In[139]:


pm.Matplot.summary_plot(expected_clearance, custom_labels=rate_labels, 
                        main='Expected clearance (%)', vline_pos=-1)


# In[123]:


μ.summary()


# In[124]:


σ.summary()