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]:
%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
Out[5]:
array(['propranolol', 'atenolol', 'Timolol maleate 0.5% gel', 'placebo',
       'triamcilone', 'prednisolone', 'control', 'control - no treatment',
       'timolol', 'Imiquimod', 'bleomycin', 'prednisone',
       'corticosteroids', 'No treatment - observation', 'triamcinolone',
       'nadolol', 'Methylprednisolone'], dtype=object)
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
Out[9]:
{'Control (not treated with propranolol)': 'control',
 'Imiquimod': 'imiquimod',
 'Methylprednisolone': 'methylprednisolone',
 'No treatment - observation': 'control',
 'Timolol maleate 0.5% gel': 'timolol',
 'atenolol': 'atenolol',
 'control': 'control',
 'control - no treatment': 'control',
 'corticosteroids': 'oral steroid',
 'methyprednisolone': 'methylprednisolone',
 'nadolol': 'nadolol',
 'nonpropranolol': 'control',
 'placebo': 'control',
 'prednisolone': 'oral steroid',
 'prednisone': 'oral steroid',
 'propranolol': 'propranolol',
 'timolol': 'timolol',
 'triamcilone': 'triamcilone',
 'triamcilone acetonide': 'triamcilone',
 'triamcinolone': 'triamcilone'}
In [10]:
kq2_outcomes['intervention'] = kq2_outcomes.Arm.replace(intervention_lookup)
kq2_outcomes['intervention'].value_counts()
Out[10]:
propranolol           34
control               15
triamcilone            9
timolol                9
oral steroid           8
atenolol               2
nadolol                1
imiquimod              1
methylprednisolone     1
dtype: int64
In [11]:
unique_interventions = kq2_outcomes.intervention.unique()
n_interventions = len(unique_interventions) - 1
n_interventions
Out[11]:
8
In [12]:
unique_interventions
Out[12]:
array(['propranolol', 'atenolol', 'timolol', 'control', 'triamcilone',
       'oral steroid', 'imiquimod', 'nadolol', 'methylprednisolone'], dtype=object)

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
Out[14]:
17
In [15]:
unique_studies
Out[15]:
array([   5,  110, 1263,   43,   13,  445, 3740,   72,  112,  309,  321,
        402,  438, 3451, 3522, 3723, 1180])

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()
Out[16]:
Dosage form    Dosage amount  intervention                 
intralesional  mg/kg          triamcilone         1 to 5        3
               mg/ml          propranolol         1            10
                              triamcilone         40            6
intravenous    mg/kg          methylprednisolone  30            1
oral           mg/kg          atenolol            1             2
                              nadolol             up to 4.0     1
                              oral steroid        2.0           5
                                                  4.0           2
                                                  2.8           1
                              propranolol         2.0          17
                                                  2 to 3        1
                                                  3             1
                                                  2.7           1
topical        BID            propranolol         0.01          4
                              timolol             0.0025        3
               Drop           control             1             1
                              timolol             1             1
               NR             timolol             NR            4
dtype: int64
In [17]:
kq2_outcomes.groupby(['intervention', 'Dosage form']).REFID.count()
Out[17]:
intervention        Dosage form  
atenolol            oral              2
control             topical           1
imiquimod           topical           1
methylprednisolone  intravenous       1
nadolol             oral              1
oral steroid        oral              8
propranolol         intralesional    10
                    oral             20
                    topical           4
timolol             topical           9
triamcilone         intralesional     9
Name: REFID, dtype: int64

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
Out[20]:
(8,)
In [21]:
σ.value.shape
Out[21]:
(9,)

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]
Out[24]:
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
0 5 propranolol 2 mg/kg oral 10 Complete response 6 (60) range NaN NaN 99 100 6
1 5 atenolol 1 mg/kg oral 13 Complete response 7 (53.8) range NaN NaN 99 100 7
2 5 propranolol 2 mg/kg oral 10 Partial response 4 (40) range NaN NaN 1 99 4
3 5 atenolol 1 mg/kg oral 13 Partial response 6 (46.1) range NaN NaN 1 99 6
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]
Out[28]:
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
4 110 Timolol maleate 0.5% gel 1 Drop topical 15 Volume reduced by ≥ 5% 60 range NaN NaN 5 100 9
5 110 placebo 1 Drop topical 18 Volume reduced by ≥ 5% 11 range NaN NaN 5 100 2
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]
Out[32]:
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
6 1263 triamcilone 1 to 5 mg/kg intralesional 25 Lesion size reduction > 50% 11 (44) range NaN NaN 50 100 11
7 1263 prednisolone 2 mg/kg oral 25 Lesion size reduction > 50% 8 (32) range NaN NaN 50 100 8
8 1263 control no treatment NaN NaN 25 Lesion size reduction > 50% 0 range NaN NaN 50 100 0
9 1263 triamcilone 1 to 5 mg/kg intralesional 25 Lesion size reduction < 50% 8 (32) range NaN NaN 10 50 8
10 1263 prednisolone 2 mg/kg oral 25 Lesion size reduction < 50% 11 (44) range NaN NaN 10 50 11
11 1263 control no treatment NaN NaN 25 Lesion size reduction < 50% 1 (4) range NaN NaN 10 50 1
12 1263 triamcilone 1 to 5 mg/kg intralesional 25 Lesion size reduction little or none 6 (24) range NaN NaN 0 10 6
13 1263 prednisolone 2 mg/kg oral 25 Lesion size reduction little or none 6 (24) range NaN NaN 0 10 6
14 1263 control no treatment NaN NaN 25 Lesion size reduction little or none 19 (76) range NaN NaN 0 10 19
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]
Out[37]:
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
15 43 propranolol 2 mg/kg oral 15 Excellent response (complete resolution) 9 (60) range NaN NaN 99 100 9
16 43 propranolol 0.01 BID topical 15 Excellent response (complete resolution) 3 (20) range NaN NaN 99 100 3
17 43 propranolol 1 mg/ml intralesional 15 Excellent response (complete resolution) 2 (13.3) range NaN NaN 99 100 2
18 43 propranolol 2 mg/kg oral 15 Good response (≥50% reduction in size of IH) 2 (13.3) range NaN NaN 50 99 2
19 43 propranolol 0.01 BID topical 15 Good response (≥50% reduction in size of IH) 5 (33.3) range NaN NaN 50 99 5
20 43 propranolol 1 mg/ml intralesional 15 Good response (≥50% reduction in size of IH) 3 (20) range NaN NaN 50 99 3
21 43 propranolol 2 mg/kg oral 15 Fair response (<50% reduction in size of IH) 1 (6.7) range NaN NaN 1 50 1
22 43 propranolol 0.01 BID topical 15 Fair response (<50% reduction in size of IH) 2 (13.3) range NaN NaN 1 50 2
23 43 propranolol 1 mg/ml intralesional 15 Fair response (<50% reduction in size of IH) 3 (20) range NaN NaN 1 50 3
24 43 propranolol 2 mg/kg oral 15 Poor response (no response, worsening of IH) 3 (20) scalar NaN NaN 0 1 3
25 43 propranolol 0.01 BID topical 15 Poor response (no response, worsening of IH) 5 (33.3) scalar NaN NaN 0 1 5
26 43 propranolol 1 mg/ml intralesional 15 Poor response (no response, worsening of IH) 7 (46.7) scalar NaN NaN 0 1 7
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]
Out[42]:
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
27 13 propranolol 2 mg/kg oral 9 Change in TSA, square millimeters 0.64 [0.44-0.83] estimate 0.64 Back calculate NaN NaN NaN
28 13 prednisolone 2 mg/kg oral 6 Change in TSA, square millimeters 0.41 [0.10-0.72] estimate 0.41 Back calculate NaN NaN NaN
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]
Out[48]:
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
29 445 propranolol 2 mg/kg oral 18 Change in volume -60 estimate 0.600 NaN NaN NaN 18
30 445 placebo NaN NaN NaN 15 Change in volume -14.1 estimate 0.147 NaN NaN NaN 15

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]
Out[49]:
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
31 3740 propranolol 3 mg/kg oral 101 Complete or nearly complete resolution 61 (60) range NaN NaN 90 100 61
32 3740 placebo NaN NaN NaN 25 Complete or nearly complete resolution 2 (4) range NaN NaN 90 100 2
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]
Out[53]:
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
33 72 propranolol 2 mg/kg oral 31 Significant involution (>50%) 28 range NaN NaN 50 100 28
34 72 propranolol 2 mg/kg oral 31 Some involution (11%-50%) 0 range NaN NaN 11 50 0
35 72 propranolol 2 mg/kg oral 31 No involution (≤10%) 3 range NaN NaN 0 10 3
36 72 control - no treatment NaN NaN NaN 14 Significant involution (>50%) 4 range NaN NaN 50 100 4
37 72 control - no treatment NaN NaN NaN 14 Some involution (11%-50%) 2 range NaN NaN 11 50 2
38 72 control - no treatment NaN NaN NaN 14 No involution (≤10%) 6 range NaN NaN 0 10 6
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]
Out[57]:
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
39 112 timolol 0.005 NaN topical 20 Visual analog scale NR IPD NaN NaN NaN NaN NR
40 112 Imiquimod 0.05 NaN topical 20 Visual analog scale NR IPD NaN NaN NaN NaN NR
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()
Out[60]:
77.5
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]
Out[62]:
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
44 309 propranolol 2 mg/kg oral 10 Complete involution (> 90% response) 2 range NaN NaN 90 100 2
45 309 propranolol 2 mg/kg oral 10 reduction in size of 75-90 % 4 range NaN NaN 75 90 4
46 309 propranolol 2 mg/kg oral 10 reduction in size of 50 to 75 % 3 range NaN NaN 50 75 3
47 309 propranolol 2 mg/kg oral 10 reduction in size of <25 % 1 range NaN NaN 0 25 1
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]
Out[65]:
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
48 321 timolol 0.0025 BID topical 13 reduction in size of >50 % 61.5 range NaN NaN 50 100 8
49 321 timolol 0.0025 BID topical 13 reduction in size of 0-50 % 30.8 range NaN NaN 1 50 4
50 321 timolol 0.0025 BID topical 13 enlarged in size NaN range NaN NaN 0 1 1
51 321 control NaN NaN NaN 10 reduction in size of >50 % 0 range NaN NaN 50 100 0
52 321 control NaN NaN NaN 10 reduction in size of 0-50 % 10 range NaN NaN 1 50 1
53 321 control NaN NaN NaN 10 enlarged in size NaN range NaN NaN 0 1 9
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]
Out[69]:
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
54 402 prednisone 2.8 mg/kg oral 12 Visual Analog Scale - mean mm 44.82 estimate 44.82 12.21 NaN NaN NaN
55 402 propranolol 2.7 mg/kg oral 12 Visual Analog Scale - mean mm 78.73 estimate 78.73 22.47 NaN NaN NaN
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]
Out[74]:
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
56 438 propranolol 2 mg/kg oral 68 Clearance ≥ 75% 56 range NaN NaN 75 100 56
57 438 propranolol 2 mg/kg oral 68 Clearance <75% 12 range NaN NaN 0 75 12
58 438 corticosteroids 4 mg/kg oral 42 Clearance ≥ 75% 12 range NaN NaN 75 100 12
59 438 corticosteroids 4 mg/kg oral 42 Clearance <75% 30 range NaN NaN 0 75 30
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]
Out[78]:
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
60 3451 timolol NR NR topical 101 Response - ineffective 7.9 range NaN NaN 0 1 8
61 3451 timolol NR NR topical 101 Response - controlled growth 35.6 range NaN NaN 0 1 36
62 3451 timolol NR NR topical 101 Response - promoted regression 56.4 range NaN NaN 1 99 57
63 3451 timolol NR NR topical 101 Complete regression 12 range NaN NaN 99 100 12
64 3451 No treatment - observation NaN NaN NaN 23 Response - ineffective 65.2 range NaN NaN 0 1 15
65 3451 No treatment - observation NaN NaN NaN 23 Response - controlled growth 30.4 range NaN NaN 0 1 7
66 3451 No treatment - observation NaN NaN NaN 23 Response - promotedregression 4.3 range NaN NaN 1 99 1
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]
Out[82]:
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
67 3522 propranolol 1 mg/ml intralesional 12 Size of hemangioma (cm2) 3.6 ± 2.6 NaN NaN NaN NaN NaN NaN
68 3522 propranolol 1 mg/ml intralesional 12 Response rate - Complete resolution achieved 42 scalar NaN NaN 99 100 5
69 3522 propranolol 1 mg/ml intralesional 12 Response rate - Sustained plateau, with >50% r... 25 range NaN NaN 50 99 3
70 3522 propranolol 1 mg/ml intralesional 12 Response rate - <50% reduction in size of hema... 17 range NaN NaN 0 50 2
71 3522 propranolol 1 mg/ml intralesional 12 Resistant to treatment 17 NaN NaN NaN NaN NaN 2
72 3522 propranolol 1 mg/ml intralesional 12 Rebound growth 4 NaN NaN NaN NaN NaN 4
73 3522 triamcinolone 40 mg/ml intralesional 10 Size of hemangioma (cm2) 3.7 ± 2.5 NaN NaN NaN NaN NaN NaN
74 3522 triamcinolone 40 mg/ml intralesional 10 Response rate - Complete resolution achieved 40 range NaN NaN 99 100 4
75 3522 triamcinolone 40 mg/ml intralesional 10 Response rate - Sustained plateau, with >50% r... 20 range NaN NaN 50 99 2
76 3522 triamcinolone 40 mg/ml intralesional 10 Response rate - <50% reduction in size of hema... 20 range NaN NaN 0 50 2
77 3522 triamcinolone 40 mg/ml intralesional 10 Resistant to treatment 20 range NaN NaN 0 0 2
78 3522 triamcinolone 40 mg/ml intralesional 10 Rebound growth 3 NaN NaN NaN NaN NaN 3
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]
Out[86]:
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
79 3723 nadolol up to 4.0 mg/kg oral 10 percentage IH shrinkage 97 ± 3.05 estimate 97 3.05 NaN NaN NaN
80 3723 propranolol 2 to 3 mg/kg oral 9 percentage IH shrinkage 86 ± 14.82 estimate 86 14.82 NaN NaN NaN
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]
Out[91]:
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
81 1180 Methylprednisolone 30 mg/kg intravenous 10 Change in size measured using visual acuity sc... -1.5 (-35 to 22) difference -1.5 NaN -35 22 NaN
82 1180 prednisolone 2 mg/kg oral 10 Change in size measured using VAS at AGE 1 year 50 (35 to 67) difference 50.0 NaN 35 67 NaN
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)
 [-                 2%                  ] 5282 of 200000 complete in 104.1 sec
In [113]:
M.sample(iterations, burn)
In [114]:
best.summary()
best:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	0.0              0.0              0.0                    [ 0.  0.]
	0.001            0.037            0.001                  [ 0.  0.]
	0.134            0.341            0.028                  [ 0.  1.]
	0.137            0.344            0.03                   [ 0.  1.]
	0.0              0.0              0.0                    [ 0.  0.]
	0.0              0.0              0.0                    [ 0.  0.]
	0.0              0.0              0.0                    [ 0.  0.]
	0.0              0.0              0.0                    [ 0.  0.]
	0.727            0.445            0.039                  [ 0.  1.]
	0.0              0.0              0.0                    [ 0.  0.]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             0.0            0.0           1.0
	0.0              0.0             0.0            0.0           1.0
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             0.0            0.0           0.0
	0.0              0.0             1.0            1.0           1.0
	0.0              0.0             0.0            0.0           0.0
	
In [125]:
pm.Matplot.plot(θ)
Plotting θ_0
Plotting θ_1
Plotting θ_2
Plotting θ_3
Plotting θ_4
Plotting θ_5
Plotting θ_6
Plotting θ_7
In [115]:
pm.Matplot.plot(σ)
Plotting σ_0
Plotting σ_1
Plotting σ_2
Plotting σ_3
Plotting σ_4
Plotting σ_5
Plotting σ_6
Plotting σ_7
Plotting σ_8
In [116]:
μ.summary()
μ:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	-2.022           0.682            0.067            [-3.296 -1.043]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	-3.321           -2.571          -1.825         -1.478        -1.048
	
In [117]:
θ.summary()
θ:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	5.31             0.772            0.071            [ 3.938  6.841]
	5.069            1.371            0.133            [ 2.412  7.345]
	3.778            0.486            0.045            [ 2.798  4.732]
	3.154            0.441            0.037            [ 2.326  3.942]
	0.847            0.892            0.087            [-0.587  2.402]
	2.679            0.608            0.056            [ 1.503  3.942]
	6.607            1.275            0.124            [ 4.273  8.866]
	1.195            0.976            0.096            [-0.144  3.149]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	3.821            4.795           5.26           5.875         6.763
	2.075            4.049           5.261          6.226         7.076
	2.852            3.444           3.779          4.115         4.82
	2.376            2.844           3.133          3.449         4.078
	-0.696           0.087           0.818          1.622         2.358
	1.457            2.261           2.718          3.058         3.916
	4.38             5.689           6.525          7.469         9.13
	-0.135           0.359           1.073          1.881         3.264
	
In [118]:
pl.hist(pm.invlogit(μ.trace() + θ.trace()[:, 5]))
Out[118]:
(array([   53.,   199.,   607.,  1059.,  1206.,  1898.,  1904.,  1560.,
         1141.,   373.]),
 array([ 0.2064767 ,  0.28074814,  0.35501958,  0.42929102,  0.50356245,
         0.57783389,  0.65210533,  0.72637677,  0.80064821,  0.87491964,
         0.94919108]),
 <a list of 10 Patch objects>)
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')
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
In [138]:
pm.Matplot.summary_plot(θ, custom_labels=sigma_labels, 
                        vline_pos=-10, main='Expected effect size (relative to control)')
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
In [121]:
pm.Matplot.summary_plot(predicted_clearance, custom_labels=rate_labels, vline_pos=-1)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
In [139]:
pm.Matplot.summary_plot(expected_clearance, custom_labels=rate_labels, 
                        main='Expected clearance (%)', vline_pos=-1)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
In [123]:
μ.summary()
μ:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	-2.022           0.682            0.067            [-3.296 -1.043]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	-3.321           -2.571          -1.825         -1.478        -1.048
	
In [124]:
σ.summary()
σ:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	2.651            0.182            0.014            [ 2.282  3.007]
	1.842            0.605            0.05             [ 0.431  2.776]
	1.503            0.033            0.003            [ 1.442  1.574]
	2.364            0.314            0.026            [ 1.852  2.951]
	1.363            0.15             0.012            [ 1.102  1.687]
	0.757            0.131            0.011            [ 0.515  0.991]
	1.269            0.989            0.079            [ 0.01   3.137]
	0.375            0.339            0.032            [ 0.04   1.057]
	1.106            0.122            0.011            [ 0.885  1.316]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	2.32             2.542           2.635          2.767         3.056
	0.536            1.422           1.909          2.212         2.944
	1.443            1.482           1.501          1.525         1.575
	1.908            2.148           2.31           2.56          3.035
	1.129            1.251           1.344          1.45          1.724
	0.523            0.664           0.747          0.859         1.018
	0.061            0.442           1.076          1.875         3.683
	0.042            0.106           0.272          0.52          1.248
	0.892            1.012           1.097          1.186         1.349