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.
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.
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})$$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.
%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
sb.set_style("white")
# Set random number seed
np.random.seed(42)
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.
unique_outcomes = kq2_outcomes['Outcome'].unique()
unique_interventions = kq2_outcomes['Arm'].unique()
unique_interventions
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)
intervention_categories = pd.read_csv('interventions.csv', index_col=0)
Remove unusable interventions
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
intervention_lookup = intervention_categories[intervention_categories.Label!='DROP'].to_dict()['Label']
intervention_lookup
{'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'}
kq2_outcomes['intervention'] = kq2_outcomes.Arm.replace(intervention_lookup)
kq2_outcomes['intervention'].value_counts()
propranolol 34 control 15 triamcilone 9 timolol 9 oral steroid 8 atenolol 2 nadolol 1 imiquimod 1 methylprednisolone 1 dtype: int64
unique_interventions = kq2_outcomes.intervention.unique()
n_interventions = len(unique_interventions) - 1
n_interventions
8
unique_interventions
array(['propranolol', 'atenolol', 'timolol', 'control', 'triamcilone', 'oral steroid', 'imiquimod', 'nadolol', 'methylprednisolone'], dtype=object)
Define indices to intervention elements
PROPRANOLOL = 0
ATENOLOL = 1
TIMOLOL = 2
TRIAMCILONE = 3
ORAL_STEROID = 4
IMIQUIMOD = 5
NADOLOL = 6
METHYLPREDNISOLONE = 7
Number of unique studies
unique_studies = kq2_outcomes.REFID.unique()
n_studies = len(unique_studies)
n_studies
17
unique_studies
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.
kq2_outcomes.groupby(['Dosage form','Dosage amount','intervention'])['Dosage'].value_counts()
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
kq2_outcomes.groupby(['intervention', 'Dosage form']).REFID.count()
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
Study random effect. Might want to use a robust distribution, like Cauchy or t.
σ_ϵ = pm.Uniform('σ_ϵ', 0, 1000, value=10)
τ_ϵ = σ_ϵ ** -2
ϵ = pm.Normal('ϵ', 0, τ_ϵ, value=np.zeros(n_studies))
Parameters for each intervention
μ = 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
θ.value.shape
(8,)
σ.value.shape
(9,)
Alternate mode of delivery effects for propanolol (oral is baseline)
ϕ = pm.Normal('ϕ', 0, 0.001, value=0)
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']
kq2_outcomes[kq2_outcomes.REFID==5][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid5 = np.where(unique_studies==5)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==110][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid110 = np.where(unique_studies==110)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==1263][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid1263 = np.where(unique_studies==1263)[0][0]
@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)
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==43][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid43 = np.where(unique_studies==43)[0][0]
@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)
# @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)
@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)
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.
kq2_outcomes[kq2_outcomes.REFID==13][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid13 = np.where(unique_studies==13)[0][0]
_propranolol_sd = 9 * (0.83 - 0.44) / 3.92
_prednisolone_sd = 6 * (0.72 - 0.10) / 3.92
propranolol_13_var = (_propranolol_sd**2) / (0.64 * (1-0.64))**2
prednisolone_13_var = (_prednisolone_sd**2) / (0.41 * (1-0.41))**2
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==445][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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.
kq2_outcomes[kq2_outcomes.REFID==3740][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid3740 = np.where(unique_studies==3740)[0][0]
@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)
@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)
The paper claims n=14
for controls, but only reports 12 outcomes.
kq2_outcomes[kq2_outcomes.REFID==72][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid72 = np.where(unique_studies==72)[0][0]
@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)
@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 et al. report individual patient data, which we incorporate directly into the meta-analysis.
kq2_outcomes[kq2_outcomes.REFID==112][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid112 = np.where(unique_studies==112)[0][0]
obs_imiquimod_112 = kq2_ipd[kq2_ipd.Arm=='imiquimod'].Outcome.values/100.
obs_timolol_112 = kq2_ipd[kq2_ipd.Arm=='timolol'].Outcome.values/100.
kq2_ipd[kq2_ipd.Arm=='imiquimod'].Outcome.median()
77.5
@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)
kq2_outcomes[kq2_outcomes.REFID==309][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid309 = np.where(unique_studies==309)[0][0]
@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)
kq2_outcomes[kq2_outcomes.REFID==321][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid321 = np.where(unique_studies==321)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==402][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid402 = np.where(unique_studies==402)[0][0]
Find SD on the logit scale using the delta method
propanolol_402_var = (0.2247**2) / (0.7873 * (1-0.7873))**2
prednisone_402_var = (0.1221**2) / (0.4482 * (1-0.4482))**2
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==438][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid438 = np.where(unique_studies==438)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==3451][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid3451 = np.where(unique_studies==3451)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==3522][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid3522 = np.where(unique_studies==3522)[0][0]
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==3723][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid3723 = np.where(unique_studies==3723)[0][0]
propanolol_3723_var = (0.1482**2) / (0.86 * (1-0.86))**2
nadolol_3723_var = (0.0305**2) / (0.97 * (1-0.97))**2
@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)
@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)
kq2_outcomes[kq2_outcomes.REFID==1180][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 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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 |
refid1180 = np.where(unique_studies==1180)[0][0]
Assuming the width of the interquartile range is 1.35 SD.
prednisolone_sd_1180 = (0.67 - 0.35) / 1.35
methylprednisolone_sd_1180 = (0.22 + 0.35) / 1.35
methylprednisolone_1180_latent = pm.Normal('methylprednisolone_1180_latent', -0.015,
methylprednisolone_sd_1180**-2, size=10)
methylprednisolone_1180_trunc = pm.Lambda('methylprednisolone_1180_trunc',
lambda y=methylprednisolone_1180_latent: np.clip(y, 1e-6, 1-1e-6))
prednisolone_1180_latent = pm.Normal('prednisolone_1180_latent', 0.5,
prednisolone_sd_1180**-2, size=10)
prednisolone_1180_trunc = pm.Lambda('prednisolone_1180_trunc',
lambda y=prednisolone_1180_latent: np.clip(y, 1e-6, 1-1e-6)).value
@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])
@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
@pm.deterministic
def expected_clearance(μ=μ, θ=θ, ϕ=ϕ):
p = np.append(pm.invlogit(μ), pm.invlogit(μ + θ[0] + ϕ))
return np.append(p, pm.invlogit(μ + θ))
@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)
best = pm.Lambda('best', lambda b=expected_clearance: (b==b.max()).astype(int))
rate_labels = ['Control',
'Intralesional propanolol',
'Oral propanolol',
'Atenolol',
'Timolol',
'Triamcilone',
'Oral steroid',
'Imiquimod',
'Nadolol',
'Methylprednisolone']
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
M.sample(iterations, burn)
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
pm.Matplot.plot(θ)
Plotting θ_0 Plotting θ_1 Plotting θ_2 Plotting θ_3 Plotting θ_4 Plotting θ_5 Plotting θ_6 Plotting θ_7
pm.Matplot.plot(σ)
Plotting σ_0 Plotting σ_1 Plotting σ_2 Plotting σ_3 Plotting σ_4 Plotting σ_5 Plotting σ_6 Plotting σ_7 Plotting σ_8
μ.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
θ.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
pl.hist(pm.invlogit(μ.trace() + θ.trace()[:, 5]))
(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>)
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.
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.
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.
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.
μ.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
σ.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