metafopy showcase¶

This notebook shows what meta-analysis looks like when it lives in Python and Jupyter. Every number comes from metafopy; the charts are Altair (Vega-Lite) charts that stay interactive on the static docs site, because the fits are computed here in Python and the browser only filters them. The statistics are those of the R package metafor (Viechtbauer, 2010, doi:10.18637/jss.v036.i03); the plots are free to look different, see Plots compared with metafor for the matplotlib versions that mirror R.

The story: a main thread on the BCG vaccine trials (Colditz et al., 1994, doi:10.1001/jama.1994.03520050036011), then two side trips: a multilevel model for schools nested in districts (Konstantopoulos, 2011, doi:10.1002/jrsm.35) and publication bias in passive smoking and lung cancer (Hackshaw et al., 1997, doi:10.1136/bmj.315.7114.980), and a generated write-up with exported figures.

Hover for tooltips, click a study to highlight it in its plot, drag in the Baujat plot to brush, and use the controls above the charts.

In [1]:

import pathlib
import tempfile

import altair
import IPython.display
import numpy
import pandas

import metafopy
import metafopy.interactive
from metafopy import cumul, datasets, escalc, fsn, gosh, influence, predict, ranef, regtest, rma, rma_mv, selmodel, trimfill

GRAPHITE = metafopy.interactive.forest.TEXT_COLOUR
BLUE = metafopy.interactive.forest.STUDY_COLOUR
GREEN = '#2AA198'
YELLOW = '#B58900'
RED = metafopy.interactive.forest.POOLED_COLOUR
MUTED = metafopy.interactive.forest.REFLINE_COLOUR

altair.renderers.enable('html', embed_options={'actions': False})
altair.data_transformers.disable_max_rows()
CONTROLS_ON_TOP = """
<style>
div.vega-embed { display: flex; flex-direction: column-reverse; align-items: center; width: 100%; }
div.vega-embed .vega-bindings { display: flex; flex-wrap: wrap; justify-content: center; align-items: center; gap: 6px 28px;
  padding: 8px 12px; margin-bottom: 10px; border: 1px solid #DDE3E8; border-radius: 6px; font-size: 14px; }
div.vega-embed .vega-bind { display: flex; align-items: center; gap: 8px; }
div.vega-embed .vega-bind label { display: inline-flex; align-items: center; gap: 8px; }
div.vega-embed .vega-bind-name { font-weight: 600; }
</style>
"""
IPython.display.display(IPython.display.HTML(CONTROLS_ON_TOP))


def bare(chart):
    """A library chart without its top-level configuration, so it can sit in a concatenation that carries one."""
    out = chart.copy(deep=True)
    out.config = altair.Undefined
    return out


print(f'metafopy {metafopy.__version__}, altair {altair.__version__}')

import pathlib import tempfile import altair import IPython.display import numpy import pandas import metafopy import metafopy.interactive from metafopy import cumul, datasets, escalc, fsn, gosh, influence, predict, ranef, regtest, rma, rma_mv, selmodel, trimfill GRAPHITE = metafopy.interactive.forest.TEXT_COLOUR BLUE = metafopy.interactive.forest.STUDY_COLOUR GREEN = '#2AA198' YELLOW = '#B58900' RED = metafopy.interactive.forest.POOLED_COLOUR MUTED = metafopy.interactive.forest.REFLINE_COLOUR altair.renderers.enable('html', embed_options={'actions': False}) altair.data_transformers.disable_max_rows() CONTROLS_ON_TOP = """ <style> div.vega-embed { display: flex; flex-direction: column-reverse; align-items: center; width: 100%; } div.vega-embed .vega-bindings { display: flex; flex-wrap: wrap; justify-content: center; align-items: center; gap: 6px 28px; padding: 8px 12px; margin-bottom: 10px; border: 1px solid #DDE3E8; border-radius: 6px; font-size: 14px; } div.vega-embed .vega-bind { display: flex; align-items: center; gap: 8px; } div.vega-embed .vega-bind label { display: inline-flex; align-items: center; gap: 8px; } div.vega-embed .vega-bind-name { font-weight: 600; } </style> """ IPython.display.display(IPython.display.HTML(CONTROLS_ON_TOP)) def bare(chart): """A library chart without its top-level configuration, so it can sit in a concatenation that carries one.""" out = chart.copy(deep=True) out.config = altair.Undefined return out print(f'metafopy {metafopy.__version__}, altair {altair.__version__}')

metafopy 0.0.1, altair 6.3.0

1. The data and a first model¶

escalc turns the 2x2 tables of the 13 trials into log risk ratios and sampling variances; rma fits the random-effects model. Results render as tables here (the printed text of R's summary() is still what print(fit) gives).

In [2]:

bcg = datasets.load('dat.bcg')
es = escalc('RR', ai='tpos', bi='tneg', ci='cpos', di='cneg', data=bcg)
es['study'] = bcg['author'] + ', ' + bcg['year'].astype(str)
es[['study', 'year', 'ablat', 'alloc', 'yi', 'vi']].head()

bcg = datasets.load('dat.bcg') es = escalc('RR', ai='tpos', bi='tneg', ci='cpos', di='cneg', data=bcg) es['study'] = bcg['author'] + ', ' + bcg['year'].astype(str) es[['study', 'year', 'ablat', 'alloc', 'yi', 'vi']].head()

Out[2]:

	study	year	ablat	alloc	yi	vi
0	Aronson, 1948	1948	44	random	-0.889311	0.325585
1	Ferguson & Simes, 1949	1949	55	random	-1.585389	0.194581
2	Rosenthal et al, 1960	1960	42	random	-1.348073	0.415368
3	Hart & Sutherland, 1977	1977	52	random	-1.441551	0.020010
4	Frimodt-Moller et al, 1973	1973	13	alternate	-0.217547	0.051210

In [3]:

fit = rma(yi='yi', vi='vi', data=es)
fit

fit = rma(yi='yi', vi='vi', data=es) fit

Out[3]:

Random-Effects Model (k = 13; tau^2 estimator: REML)

	estimate	se	zval	pval	ci.lb	ci.ub
intrcpt	-0.7145	0.1798	-3.974	7.054e-05	-1.067	-0.3622

k = 13; τ² = 0.3132; I² = 92.22%; H² = 12.86; Q(df = 12) = 152.23, p-val <.0001

The pooled estimate is on the log scale. A prediction interval says where the effect of a new trial would fall, which is much wider than the confidence interval when the trials disagree as much as these do.

In [4]:

pred = predict(fit, transf=numpy.exp)
print(
    f'Pooled risk ratio {float(pred.pred[0]):.2f}, 95% CI [{float(pred.ci_lb[0]):.2f}, {float(pred.ci_ub[0]):.2f}], '
    f'95% prediction interval [{float(pred.pi_lb[0]):.2f}, {float(pred.pi_ub[0]):.2f}]'
)

pred = predict(fit, transf=numpy.exp) print( f'Pooled risk ratio {float(pred.pred[0]):.2f}, 95% CI [{float(pred.ci_lb[0]):.2f}, {float(pred.ci_ub[0]):.2f}], ' f'95% prediction interval [{float(pred.pi_lb[0]):.2f}, {float(pred.pi_ub[0]):.2f}]' )

Pooled risk ratio 0.49, 95% CI [0.34, 0.70], 95% prediction interval [0.15, 1.55]

2. Everything a study contributes, in one frame¶

The charts below share one tidy frame with a row per study: effect and interval, weight, the coordinates of the funnel and Baujat plots, and the influence diagnostics of influence().

In [5]:

yi = numpy.asarray(es['yi'], dtype=float)
vi = numpy.asarray(es['vi'], dtype=float)
inf = influence(fit).to_frame()

fixed = rma(yi='yi', vi='vi', data=es, method='EE')
theta_fixed = float(fixed.beta[0])
baujat_y = []
for i in range(len(yi)):
    keep = numpy.arange(len(yi)) != i
    loo = rma(yi=yi[keep], vi=vi[keep], method='EE')
    baujat_y.append((theta_fixed - float(loo.beta[0])) ** 2 / float(loo.se[0]) ** 2)

weight = 1 / (vi + fit.tau2)
studies = pandas.DataFrame(
    {
        'study': es['study'],
        'alloc': bcg['alloc'],
        'ablat': bcg['ablat'],
        'n': bcg[['tpos', 'tneg', 'cpos', 'cneg']].sum(axis=1),
        'rr': numpy.exp(yi),
        'lo': numpy.exp(yi - 1.96 * numpy.sqrt(vi)),
        'hi': numpy.exp(yi + 1.96 * numpy.sqrt(vi)),
        'yi': yi,
        'sei': numpy.sqrt(vi),
        'weight': 100 * weight / weight.sum(),
        'baujat_x': (yi - theta_fixed) ** 2 / vi,
        'baujat_y': baujat_y,
        'cook_d': numpy.asarray(inf['cook.d'], dtype=float),
        'rstudent': numpy.asarray(inf['rstudent'], dtype=float),
    }
)
order = list(studies.sort_values('ablat')['study'])
studies.round(3).head()

yi = numpy.asarray(es['yi'], dtype=float) vi = numpy.asarray(es['vi'], dtype=float) inf = influence(fit).to_frame() fixed = rma(yi='yi', vi='vi', data=es, method='EE') theta_fixed = float(fixed.beta[0]) baujat_y = [] for i in range(len(yi)): keep = numpy.arange(len(yi)) != i loo = rma(yi=yi[keep], vi=vi[keep], method='EE') baujat_y.append((theta_fixed - float(loo.beta[0])) ** 2 / float(loo.se[0]) ** 2) weight = 1 / (vi + fit.tau2) studies = pandas.DataFrame( { 'study': es['study'], 'alloc': bcg['alloc'], 'ablat': bcg['ablat'], 'n': bcg[['tpos', 'tneg', 'cpos', 'cneg']].sum(axis=1), 'rr': numpy.exp(yi), 'lo': numpy.exp(yi - 1.96 * numpy.sqrt(vi)), 'hi': numpy.exp(yi + 1.96 * numpy.sqrt(vi)), 'yi': yi, 'sei': numpy.sqrt(vi), 'weight': 100 * weight / weight.sum(), 'baujat_x': (yi - theta_fixed) ** 2 / vi, 'baujat_y': baujat_y, 'cook_d': numpy.asarray(inf['cook.d'], dtype=float), 'rstudent': numpy.asarray(inf['rstudent'], dtype=float), } ) order = list(studies.sort_values('ablat')['study']) studies.round(3).head()

Out[5]:

	study	alloc	ablat	n	rr	lo	hi	yi	sei	weight	baujat_x	baujat_y	cook_d	rstudent
0	Aronson, 1948	random	44	262	0.411	0.134	1.257	-0.889	0.571	5.059	0.647	0.003	0.002	-0.218
1	Ferguson & Simes, 1949	random	55	609	0.205	0.086	0.486	-1.585	0.441	6.365	6.857	0.058	0.113	-1.292
2	Rosenthal et al, 1960	random	42	451	0.260	0.073	0.919	-1.348	0.644	4.436	2.028	0.008	0.026	-0.755
3	Hart & Sutherland, 1977	random	52	26465	0.237	0.179	0.312	-1.442	0.141	9.699	51.107	4.563	0.229	-1.451
4	Frimodt-Moller et al, 1973	alternate	13	10877	0.804	0.516	1.254	-0.218	0.226	8.868	0.884	0.029	0.076	0.848

3. Fit every estimator, with and without each trial¶

The pooled estimate depends on two choices: which estimator of the between-trial variance tau^2 you use (metafopy has seven, see Veroniki et al. (2016)), and whether you leave one trial out to see if it drives the result. Here every combination is fitted once, in Python: 7 estimators x 14 options (leave out none, or one of the 13 trials) = 98 fits. Each row of the pooled frame is one fit, with its pooled risk ratio, intervals, tau^2 and I^2. The next cell prints the range of the pooled risk ratio (est) and of tau^2 over the 14 options for each estimator; section 4 turns the same frame into controls.

In [6]:

estimators = ['DL', 'EB', 'HE', 'ML', 'PM', 'REML', 'SJ']
rows = []
for method in estimators:
    for drop in ['(none)', *studies['study']]:
        keep = (studies['study'] != drop).to_numpy()
        f = rma(yi=yi[keep], vi=vi[keep], method=method)
        p = predict(f)
        rows.append(
            {
                'method': method,
                'drop': drop,
                'est': float(numpy.exp(f.beta[0])),
                'lo': float(numpy.exp(f.ci_lb[0])),
                'hi': float(numpy.exp(f.ci_ub[0])),
                'pi_lo': float(numpy.exp(numpy.ravel(p.pi_lb)[0])),
                'pi_hi': float(numpy.exp(numpy.ravel(p.pi_ub)[0])),
                'tau2': float(f.tau2),
                'I2': float(f.I2),
            }
        )
pooled = pandas.DataFrame(rows)
pooled.groupby('method')[['est', 'tau2']].agg(['min', 'max']).round(3)

estimators = ['DL', 'EB', 'HE', 'ML', 'PM', 'REML', 'SJ'] rows = [] for method in estimators: for drop in ['(none)', *studies['study']]: keep = (studies['study'] != drop).to_numpy() f = rma(yi=yi[keep], vi=vi[keep], method=method) p = predict(f) rows.append( { 'method': method, 'drop': drop, 'est': float(numpy.exp(f.beta[0])), 'lo': float(numpy.exp(f.ci_lb[0])), 'hi': float(numpy.exp(f.ci_ub[0])), 'pi_lo': float(numpy.exp(numpy.ravel(p.pi_lb)[0])), 'pi_hi': float(numpy.exp(numpy.ravel(p.pi_ub)[0])), 'tau2': float(f.tau2), 'I2': float(f.I2), } ) pooled = pandas.DataFrame(rows) pooled.groupby('method')[['est', 'tau2']].agg(['min', 'max']).round(3)

Out[6]:

	est		tau2
	min	max	min	max
method
DL	0.454	0.539	0.211	0.388
EB	0.451	0.531	0.284	0.357
HE	0.451	0.529	0.265	0.385
ML	0.453	0.537	0.226	0.320
PM	0.451	0.531	0.284	0.357
REML	0.452	0.533	0.263	0.360
SJ	0.451	0.528	0.306	0.383

4. Forest, funnel, Baujat and meta-regression plots¶

The forest and funnel plots are metafopy.interactive.forest_chart() and funnel_chart(): each takes the fit and returns an Altair chart with tooltips and click-to-highlight. The rest is composed here around them. The radio buttons and the "Leave out" drop-down above the charts pick a row of the pooled frame of section 3, the tau^2 estimator and the omitted trial; the numbers and the "Selected fit" row below them follow. Click a study (shift-click for several) in the Baujat plot to highlight it in the meta-regression, or brush a region of the Baujat plot to light up the studies in it. The bubble plot is the meta-regression on absolute latitude, the moderator that explains most of the heterogeneity (predict(..., newmods=) draws the band).

In [7]:

WIDTH = 560  # plot width that, with the axes, fills the docs page column; Vega-Lite cannot size linked plots to their container

mr = rma(yi='yi', vi='vi', mods='~ ablat', data=es)
grid = numpy.linspace(10, 60, 51)
pr = predict(mr, newmods=grid)
band = pandas.DataFrame(
    {'ablat': grid, 'pred': numpy.exp(numpy.ravel(pr.pred)), 'lo': numpy.exp(numpy.ravel(pr.ci_lb)), 'hi': numpy.exp(numpy.ravel(pr.ci_ub))}
)

by_latitude = es.sort_values('ablat')
fit_by_latitude = rma(yi='yi', vi='vi', data=by_latitude, slab=by_latitude['study'])
forest = metafopy.interactive.forest_chart(fit_by_latitude, transf=numpy.exp, refline=1.0, xlab='Risk ratio', width=WIDTH).properties(
    title='Forest plot (sorted by latitude)'
)
funnel = metafopy.interactive.funnel_chart(fit, slab=es['study'], width=WIDTH).properties(title='Funnel plot')

pick = altair.selection_point(fields=['study'], on='click', clear='dblclick', toggle='true')
brush = altair.selection_interval(encodings=['x', 'y'], name='brush')
method = altair.param(name='method', value='REML', bind=altair.binding_radio(options=estimators, name='tau² estimator  '))
drop = altair.param(name='drop', value='(none)', bind=altair.binding_select(options=['(none)', *order], name='Leave out  '))
dim = altair.condition(pick | brush, altair.value(1), altair.value(0.2))

tooltip = [
    altair.Tooltip('study:N', title='Study'),
    altair.Tooltip('rr:Q', title='Risk ratio', format='.2f'),
    altair.Tooltip('lo:Q', title='95% CI low', format='.2f'),
    altair.Tooltip('hi:Q', title='95% CI high', format='.2f'),
    altair.Tooltip('weight:Q', title='Weight %', format='.1f'),
    altair.Tooltip('ablat:Q', title='Abs. latitude'),
    altair.Tooltip('alloc:N', title='Allocation'),
    altair.Tooltip('cook_d:Q', title="Cook's distance", format='.3f'),
]
alloc_color = altair.Color('alloc:N', title='Allocation', scale=altair.Scale(domain=['random', 'alternate', 'systematic'], range=[BLUE, GREEN, YELLOW]))
base = altair.Chart(studies)

pooled_sel = altair.Chart(pooled).transform_filter((altair.datum.method == method) & (altair.datum.drop == drop))
sel_axis = altair.Axis(values=[0.05, 0.1, 0.25, 0.5, 1, 2, 4], format='~g')
sel_x = altair.X('pi_lo:Q', scale=altair.Scale(type='log', zero=False), title='Risk ratio (log scale)', axis=sel_axis)
selected_fit = (
    (
        pooled_sel.mark_rule(color=RED, strokeWidth=1).encode(x=sel_x, x2='pi_hi:Q')
        + pooled_sel.mark_rule(color=RED, strokeWidth=3).encode(x='lo:Q', x2='hi:Q')
        + pooled_sel.mark_point(shape='diamond', filled=True, size=160, color=RED).encode(x='est:Q')
    )
    .properties(width=WIDTH, height=30, title='Selected fit: estimate, 95% CI (thick) and prediction interval (thin)')
    .add_params(method, drop)
)
pooled_text = (
    pooled_sel.transform_calculate(
        txt='"RR " + format(datum.est, ".2f") + " [" + format(datum.lo, ".2f") + ", " + format(datum.hi, ".2f") + "]   PI ["'
        ' + format(datum.pi_lo, ".2f") + ", " + format(datum.pi_hi, ".2f") + "]   τ² " + format(datum.tau2, ".3f") + "   I² " + format(datum.I2, ".1f") + "%"'
    )
    .mark_text(align='left', fontSize=11, color=GRAPHITE, fontWeight='bold')
    .encode(text='txt:N')
    .properties(width=WIDTH, height=18)
)

baujat_chart = (
    base.mark_circle(size=90, stroke='white', strokeWidth=1)
    .encode(
        x=altair.X('baujat_x:Q', title='Contribution to Q'),
        y=altair.Y('baujat_y:Q', title='Influence on pooled estimate'),
        color=alloc_color,
        opacity=dim,
        tooltip=tooltip,
    )
    .add_params(pick, brush)
    .properties(width=WIDTH, height=260, title='Baujat plot (drag to brush)')
)

lat_axis = altair.Axis(values=[0.1, 0.25, 0.5, 1, 2], format='~g')
bubble_chart = (
    altair.Chart(band)
    .mark_area(color=BLUE, opacity=0.12)
    .encode(
        x=altair.X('ablat:Q', title='Absolute latitude (°)'),
        y=altair.Y('lo:Q', scale=altair.Scale(type='log', domain=[0.1, 2]), title='Risk ratio', axis=lat_axis),
        y2='hi:Q',
    )
    + altair.Chart(band).mark_line(color=BLUE, strokeWidth=2).encode(x='ablat:Q', y='pred:Q')
    + base.mark_circle(stroke='white', strokeWidth=1).encode(
        x='ablat:Q',
        y=altair.Y('rr:Q', scale=altair.Scale(type='log')),
        size=altair.Size('weight:Q', scale=altair.Scale(range=[30, 500]), legend=None),
        color=alloc_color,
        opacity=dim,
        tooltip=tooltip,
    )
).properties(width=WIDTH, height=260, title='Meta-regression on latitude (REML)')

dashboard = (
    altair.vconcat(pooled_text, selected_fit, bare(forest), bare(funnel), baujat_chart, bubble_chart, spacing=24)
    .resolve_scale(size='independent')
    .configure(font='Source Sans 3, Helvetica, Arial, sans-serif')
    .configure_view(stroke=None)
    .configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5', domainColor=MUTED, tickColor=MUTED)
    .configure_title(color=GRAPHITE, fontSize=13, anchor='start', fontWeight=600)
    .configure_legend(labelColor=GRAPHITE, titleColor=GRAPHITE, orient='top')
)
dashboard

WIDTH = 560 # plot width that, with the axes, fills the docs page column; Vega-Lite cannot size linked plots to their container mr = rma(yi='yi', vi='vi', mods='~ ablat', data=es) grid = numpy.linspace(10, 60, 51) pr = predict(mr, newmods=grid) band = pandas.DataFrame( {'ablat': grid, 'pred': numpy.exp(numpy.ravel(pr.pred)), 'lo': numpy.exp(numpy.ravel(pr.ci_lb)), 'hi': numpy.exp(numpy.ravel(pr.ci_ub))} ) by_latitude = es.sort_values('ablat') fit_by_latitude = rma(yi='yi', vi='vi', data=by_latitude, slab=by_latitude['study']) forest = metafopy.interactive.forest_chart(fit_by_latitude, transf=numpy.exp, refline=1.0, xlab='Risk ratio', width=WIDTH).properties( title='Forest plot (sorted by latitude)' ) funnel = metafopy.interactive.funnel_chart(fit, slab=es['study'], width=WIDTH).properties(title='Funnel plot') pick = altair.selection_point(fields=['study'], on='click', clear='dblclick', toggle='true') brush = altair.selection_interval(encodings=['x', 'y'], name='brush') method = altair.param(name='method', value='REML', bind=altair.binding_radio(options=estimators, name='tau² estimator ')) drop = altair.param(name='drop', value='(none)', bind=altair.binding_select(options=['(none)', *order], name='Leave out ')) dim = altair.condition(pick | brush, altair.value(1), altair.value(0.2)) tooltip = [ altair.Tooltip('study:N', title='Study'), altair.Tooltip('rr:Q', title='Risk ratio', format='.2f'), altair.Tooltip('lo:Q', title='95% CI low', format='.2f'), altair.Tooltip('hi:Q', title='95% CI high', format='.2f'), altair.Tooltip('weight:Q', title='Weight %', format='.1f'), altair.Tooltip('ablat:Q', title='Abs. latitude'), altair.Tooltip('alloc:N', title='Allocation'), altair.Tooltip('cook_d:Q', title="Cook's distance", format='.3f'), ] alloc_color = altair.Color('alloc:N', title='Allocation', scale=altair.Scale(domain=['random', 'alternate', 'systematic'], range=[BLUE, GREEN, YELLOW])) base = altair.Chart(studies) pooled_sel = altair.Chart(pooled).transform_filter((altair.datum.method == method) & (altair.datum.drop == drop)) sel_axis = altair.Axis(values=[0.05, 0.1, 0.25, 0.5, 1, 2, 4], format='~g') sel_x = altair.X('pi_lo:Q', scale=altair.Scale(type='log', zero=False), title='Risk ratio (log scale)', axis=sel_axis) selected_fit = ( ( pooled_sel.mark_rule(color=RED, strokeWidth=1).encode(x=sel_x, x2='pi_hi:Q') + pooled_sel.mark_rule(color=RED, strokeWidth=3).encode(x='lo:Q', x2='hi:Q') + pooled_sel.mark_point(shape='diamond', filled=True, size=160, color=RED).encode(x='est:Q') ) .properties(width=WIDTH, height=30, title='Selected fit: estimate, 95% CI (thick) and prediction interval (thin)') .add_params(method, drop) ) pooled_text = ( pooled_sel.transform_calculate( txt='"RR " + format(datum.est, ".2f") + " [" + format(datum.lo, ".2f") + ", " + format(datum.hi, ".2f") + "] PI ["' ' + format(datum.pi_lo, ".2f") + ", " + format(datum.pi_hi, ".2f") + "] τ² " + format(datum.tau2, ".3f") + " I² " + format(datum.I2, ".1f") + "%"' ) .mark_text(align='left', fontSize=11, color=GRAPHITE, fontWeight='bold') .encode(text='txt:N') .properties(width=WIDTH, height=18) ) baujat_chart = ( base.mark_circle(size=90, stroke='white', strokeWidth=1) .encode( x=altair.X('baujat_x:Q', title='Contribution to Q'), y=altair.Y('baujat_y:Q', title='Influence on pooled estimate'), color=alloc_color, opacity=dim, tooltip=tooltip, ) .add_params(pick, brush) .properties(width=WIDTH, height=260, title='Baujat plot (drag to brush)') ) lat_axis = altair.Axis(values=[0.1, 0.25, 0.5, 1, 2], format='~g') bubble_chart = ( altair.Chart(band) .mark_area(color=BLUE, opacity=0.12) .encode( x=altair.X('ablat:Q', title='Absolute latitude (°)'), y=altair.Y('lo:Q', scale=altair.Scale(type='log', domain=[0.1, 2]), title='Risk ratio', axis=lat_axis), y2='hi:Q', ) + altair.Chart(band).mark_line(color=BLUE, strokeWidth=2).encode(x='ablat:Q', y='pred:Q') + base.mark_circle(stroke='white', strokeWidth=1).encode( x='ablat:Q', y=altair.Y('rr:Q', scale=altair.Scale(type='log')), size=altair.Size('weight:Q', scale=altair.Scale(range=[30, 500]), legend=None), color=alloc_color, opacity=dim, tooltip=tooltip, ) ).properties(width=WIDTH, height=260, title='Meta-regression on latitude (REML)') dashboard = ( altair.vconcat(pooled_text, selected_fit, bare(forest), bare(funnel), baujat_chart, bubble_chart, spacing=24) .resolve_scale(size='independent') .configure(font='Source Sans 3, Helvetica, Arial, sans-serif') .configure_view(stroke=None) .configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5', domainColor=MUTED, tickColor=MUTED) .configure_title(color=GRAPHITE, fontSize=13, anchor='start', fontWeight=600) .configure_legend(labelColor=GRAPHITE, titleColor=GRAPHITE, orient='top') ) dashboard

/home/runner/work/metafopy/metafopy/.venv/lib/python3.14/site-packages/IPython/core/interactiveshell.py:3823: UserWarning: Automatically deduplicated selection parameter with identical configuration. If you want independent parameters, explicitly name them differently (e.g., name='param1', name='param2'). See https://github.com/vega/altair/issues/3891
  exec(code_obj, self.user_global_ns, self.user_ns)

Out[7]:

Things to try: leave out Hart & Sutherland, 1977 and watch the pooled estimate and tau^2 change; pick DL and then SJ; brush the influential trials at the top right of the Baujat plot and see where they sit in the meta-regression.

5. Cumulative meta-analysis¶

cumul() adds the trials one at a time in the order you give, here by publication year. The chart lets you change the estimator and shows how the estimate settles as the evidence accumulates.

In [8]:

years = bcg['year'].to_numpy()
cum_rows = []
for m in ['DL', 'HE', 'PM', 'REML']:
    c = cumul(rma(yi='yi', vi='vi', data=es, method=m, slab=es['study']), order=years)
    cum_rows.append(
        pandas.DataFrame(
            {
                'method': m,
                'step': numpy.arange(1, len(c.slab) + 1),
                'study': list(c.slab),
                'year': numpy.sort(years),
                'est': numpy.exp(numpy.asarray(c.estimate)),
                'lo': numpy.exp(numpy.asarray(c.ci_lb)),
                'hi': numpy.exp(numpy.asarray(c.ci_ub)),
            }
        )
    )
cumulative = pandas.concat(cum_rows, ignore_index=True)

cum_method = altair.param(name='cum_method', value='REML', bind=altair.binding_radio(options=['DL', 'HE', 'PM', 'REML'], name='tau² estimator  '))
cum_base = altair.Chart(cumulative).transform_filter(altair.datum.method == cum_method)
cum_x = altair.X('step:O', title='Trials added (in order of publication year)', axis=altair.Axis(labelAngle=0))
cum_chart = (
    (
        cum_base.mark_area(color=BLUE, opacity=0.15).encode(
            x=cum_x,
            y=altair.Y(
                'lo:Q', scale=altair.Scale(type='log', domain=[0.1, 3]), title='Pooled risk ratio', axis=altair.Axis(values=[0.1, 0.25, 0.5, 1, 2], format='~g')
            ),
            y2='hi:Q',
        )
        + cum_base.mark_line(color=BLUE, point=altair.OverlayMarkDef(filled=True, size=60)).encode(
            x=cum_x,
            y='est:Q',
            tooltip=['study:N', 'year:Q', altair.Tooltip('est:Q', format='.2f'), altair.Tooltip('lo:Q', format='.2f'), altair.Tooltip('hi:Q', format='.2f')],
        )
        + altair.Chart(pandas.DataFrame({'y': [1]})).mark_rule(color=MUTED, strokeDash=[4, 3]).encode(y='y:Q')
    )
    .properties(width=WIDTH, height=260, title='Cumulative meta-analysis')
    .add_params(cum_method)
)
cum_chart

years = bcg['year'].to_numpy() cum_rows = [] for m in ['DL', 'HE', 'PM', 'REML']: c = cumul(rma(yi='yi', vi='vi', data=es, method=m, slab=es['study']), order=years) cum_rows.append( pandas.DataFrame( { 'method': m, 'step': numpy.arange(1, len(c.slab) + 1), 'study': list(c.slab), 'year': numpy.sort(years), 'est': numpy.exp(numpy.asarray(c.estimate)), 'lo': numpy.exp(numpy.asarray(c.ci_lb)), 'hi': numpy.exp(numpy.asarray(c.ci_ub)), } ) ) cumulative = pandas.concat(cum_rows, ignore_index=True) cum_method = altair.param(name='cum_method', value='REML', bind=altair.binding_radio(options=['DL', 'HE', 'PM', 'REML'], name='tau² estimator ')) cum_base = altair.Chart(cumulative).transform_filter(altair.datum.method == cum_method) cum_x = altair.X('step:O', title='Trials added (in order of publication year)', axis=altair.Axis(labelAngle=0)) cum_chart = ( ( cum_base.mark_area(color=BLUE, opacity=0.15).encode( x=cum_x, y=altair.Y( 'lo:Q', scale=altair.Scale(type='log', domain=[0.1, 3]), title='Pooled risk ratio', axis=altair.Axis(values=[0.1, 0.25, 0.5, 1, 2], format='~g') ), y2='hi:Q', ) + cum_base.mark_line(color=BLUE, point=altair.OverlayMarkDef(filled=True, size=60)).encode( x=cum_x, y='est:Q', tooltip=['study:N', 'year:Q', altair.Tooltip('est:Q', format='.2f'), altair.Tooltip('lo:Q', format='.2f'), altair.Tooltip('hi:Q', format='.2f')], ) + altair.Chart(pandas.DataFrame({'y': [1]})).mark_rule(color=MUTED, strokeDash=[4, 3]).encode(y='y:Q') ) .properties(width=WIDTH, height=260, title='Cumulative meta-analysis') .add_params(cum_method) ) cum_chart

Out[8]:

6. GOSH: all subsets of the studies¶

Olkin, Dahabreh and Trikalinos (2012) fit the model to every possible subset of studies (here 2^13 - 1 = 8191 fits) and plot the pooled estimate against I^2. Clusters point at a subgroup of studies that behave differently. Brush a region of the scatterplot: the histogram shows the estimates in the brush, and the bars show how often each trial is in those subsets.

In [9]:

g = gosh(fit)
gosh_frame = g.res[['estimate', 'I2']].copy()
gosh_frame['estimate'] = numpy.exp(gosh_frame['estimate']).round(4)
gosh_frame['I2'] = gosh_frame['I2'].round(2)
incl = numpy.asarray(g.incl, dtype=int)
labels = list(studies['study'])
for j, label in enumerate(labels):
    gosh_frame[label] = incl[:, j]

gosh_brush = altair.selection_interval(name='gosh_brush')
scatter = (
    altair.Chart(gosh_frame)
    .mark_circle(size=8, opacity=0.25)
    .encode(
        x=altair.X('estimate:Q', scale=altair.Scale(type='log'), title='Pooled risk ratio', axis=altair.Axis(values=[0.2, 0.3, 0.5, 0.75, 1], format='~g')),
        y=altair.Y('I2:Q', title='I² (%)'),
        color=altair.condition(gosh_brush, altair.value(BLUE), altair.value(MUTED)),
    )
    .add_params(gosh_brush)
    .properties(width=WIDTH, height=260, title='All 8191 subsets (drag to brush)')
)
hist = (
    altair.Chart(gosh_frame)
    .transform_filter(gosh_brush)
    .mark_bar(color=BLUE)
    .encode(
        x=altair.X('estimate:Q', bin=altair.Bin(maxbins=30), scale=altair.Scale(type='log'), title='Pooled risk ratio in brush'),
        y=altair.Y('count()', title='Subsets'),
    )
    .properties(width=WIDTH, height=140)
)
inclusion = (
    altair.Chart(gosh_frame)
    .transform_filter(gosh_brush)
    .transform_fold(labels, as_=['study', 'in'])
    .mark_bar(color=GREEN)
    .encode(
        y=altair.Y('study:N', sort=order, title=None),
        x=altair.X('mean(in):Q', scale=altair.Scale(domain=[0, 1]), title='Share of brushed subsets containing the trial'),
    )
    .properties(width=WIDTH, height=240)
)
altair.vconcat(scatter, hist, inclusion, spacing=24).configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5').configure_title(
    color=GRAPHITE, anchor='start'
)

g = gosh(fit) gosh_frame = g.res[['estimate', 'I2']].copy() gosh_frame['estimate'] = numpy.exp(gosh_frame['estimate']).round(4) gosh_frame['I2'] = gosh_frame['I2'].round(2) incl = numpy.asarray(g.incl, dtype=int) labels = list(studies['study']) for j, label in enumerate(labels): gosh_frame[label] = incl[:, j] gosh_brush = altair.selection_interval(name='gosh_brush') scatter = ( altair.Chart(gosh_frame) .mark_circle(size=8, opacity=0.25) .encode( x=altair.X('estimate:Q', scale=altair.Scale(type='log'), title='Pooled risk ratio', axis=altair.Axis(values=[0.2, 0.3, 0.5, 0.75, 1], format='~g')), y=altair.Y('I2:Q', title='I² (%)'), color=altair.condition(gosh_brush, altair.value(BLUE), altair.value(MUTED)), ) .add_params(gosh_brush) .properties(width=WIDTH, height=260, title='All 8191 subsets (drag to brush)') ) hist = ( altair.Chart(gosh_frame) .transform_filter(gosh_brush) .mark_bar(color=BLUE) .encode( x=altair.X('estimate:Q', bin=altair.Bin(maxbins=30), scale=altair.Scale(type='log'), title='Pooled risk ratio in brush'), y=altair.Y('count()', title='Subsets'), ) .properties(width=WIDTH, height=140) ) inclusion = ( altair.Chart(gosh_frame) .transform_filter(gosh_brush) .transform_fold(labels, as_=['study', 'in']) .mark_bar(color=GREEN) .encode( y=altair.Y('study:N', sort=order, title=None), x=altair.X('mean(in):Q', scale=altair.Scale(domain=[0, 1]), title='Share of brushed subsets containing the trial'), ) .properties(width=WIDTH, height=240) ) altair.vconcat(scatter, hist, inclusion, spacing=24).configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5').configure_title( color=GRAPHITE, anchor='start' )

Out[9]:

Side trip: schools nested in districts¶

dat.konstantopoulos2011 has 56 studies of the effect of modified school calendars on achievement, from 11 districts. Studies in one district are not independent, so rma_mv adds a random effect for district and one for school within district (Konstantopoulos, 2011).

In [10]:

kdat = datasets.load('dat.konstantopoulos2011')
mv = rma_mv(yi='yi', V='vi', random='~ 1 | district/school', data=kdat)
mv

kdat = datasets.load('dat.konstantopoulos2011') mv = rma_mv(yi='yi', V='vi', random='~ 1 | district/school', data=kdat) mv

Out[10]:

Multivariate Meta-Analysis Model (k = 56; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed           factor
sigma^2.1  0.0651  0.2551     11     no         district
sigma^2.2  0.0327  0.1809     56     no  district/school

Test for Heterogeneity:
Q(df = 55) = 578.8640, p-val < .0001

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub
  0.1847  0.0846  2.1845  0.0289  0.0190  0.3504  *

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The two variance components are about equal in size to the typical sampling variance. ranef predicts the effect of each district with its prediction interval; the chart shows these over the observed study effects.

In [11]:

districts = ranef(mv)['district'].reset_index(names='district')
districts['district'] = districts['district'].astype(str)
districts['est'] = districts['intrcpt'] + float(mv.beta[0])
districts['lo'] = districts['pi.lb'] + float(mv.beta[0])
districts['hi'] = districts['pi.ub'] + float(mv.beta[0])
kstudies = kdat.assign(district=kdat['district'].astype(str), se=numpy.sqrt(kdat['vi']))
dist_order = list(districts.sort_values('est')['district'])

district_pick = altair.selection_point(fields=['district'], on='click', clear='dblclick')
dist_y = altair.Y('district:N', sort=dist_order, title='District')
dist_bars = (
    altair.Chart(districts)
    .mark_rule(strokeWidth=3, color=BLUE)
    .encode(
        x=altair.X('lo:Q', title='Standardised mean difference'),
        x2='hi:Q',
        y=dist_y,
        opacity=altair.condition(district_pick, altair.value(1), altair.value(0.3)),
    )
)
dist_pts = (
    altair.Chart(districts)
    .mark_point(filled=True, size=90, color=GRAPHITE)
    .encode(x='est:Q', y=dist_y, tooltip=['district:N', altair.Tooltip('est:Q', title='District effect', format='.2f')])
    .add_params(district_pick)
)
study_pts = (
    altair.Chart(kstudies)
    .mark_circle(size=40, color=YELLOW)
    .encode(
        x='yi:Q',
        y=dist_y,
        opacity=altair.condition(district_pick, altair.value(0.9), altair.value(0.2)),
        tooltip=['district:N', 'school:Q', 'year:Q', 'yi:Q', altair.Tooltip('se:Q', format='.2f')],
    )
)
overall = altair.Chart(pandas.DataFrame({'x': [float(mv.beta[0])]})).mark_rule(color=RED, strokeDash=[4, 3]).encode(x='x:Q')
(overall + study_pts + dist_bars + dist_pts).properties(
    width=WIDTH, height=300, title='District effects (black, with 95% prediction interval) and the school effects within them (yellow); red = overall'
).configure_title(color=GRAPHITE, fontSize=12, anchor='start').configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5')

districts = ranef(mv)['district'].reset_index(names='district') districts['district'] = districts['district'].astype(str) districts['est'] = districts['intrcpt'] + float(mv.beta[0]) districts['lo'] = districts['pi.lb'] + float(mv.beta[0]) districts['hi'] = districts['pi.ub'] + float(mv.beta[0]) kstudies = kdat.assign(district=kdat['district'].astype(str), se=numpy.sqrt(kdat['vi'])) dist_order = list(districts.sort_values('est')['district']) district_pick = altair.selection_point(fields=['district'], on='click', clear='dblclick') dist_y = altair.Y('district:N', sort=dist_order, title='District') dist_bars = ( altair.Chart(districts) .mark_rule(strokeWidth=3, color=BLUE) .encode( x=altair.X('lo:Q', title='Standardised mean difference'), x2='hi:Q', y=dist_y, opacity=altair.condition(district_pick, altair.value(1), altair.value(0.3)), ) ) dist_pts = ( altair.Chart(districts) .mark_point(filled=True, size=90, color=GRAPHITE) .encode(x='est:Q', y=dist_y, tooltip=['district:N', altair.Tooltip('est:Q', title='District effect', format='.2f')]) .add_params(district_pick) ) study_pts = ( altair.Chart(kstudies) .mark_circle(size=40, color=YELLOW) .encode( x='yi:Q', y=dist_y, opacity=altair.condition(district_pick, altair.value(0.9), altair.value(0.2)), tooltip=['district:N', 'school:Q', 'year:Q', 'yi:Q', altair.Tooltip('se:Q', format='.2f')], ) ) overall = altair.Chart(pandas.DataFrame({'x': [float(mv.beta[0])]})).mark_rule(color=RED, strokeDash=[4, 3]).encode(x='x:Q') (overall + study_pts + dist_bars + dist_pts).properties( width=WIDTH, height=300, title='District effects (black, with 95% prediction interval) and the school effects within them (yellow); red = overall' ).configure_title(color=GRAPHITE, fontSize=12, anchor='start').configure_axis(labelColor=GRAPHITE, titleColor=GRAPHITE, gridColor='#EEE8D5')

Out[11]:

Side trip: publication bias¶

dat.hackshaw1998 pools 37 studies of lung cancer risk and passive smoking. The funnel plot below marks the studies that the trim-and-fill method (Duval and Tweedie, 2000, doi:10.1111/j.0006-341X.2000.00455.x) would add to make it symmetric; hover them to see how much each would change the picture.

In [12]:

hack = datasets.load('dat.hackshaw1998')
fit_h = rma(yi='yi', vi='vi', data=hack)
tf = trimfill(fit_h)
funnel_h = (
    metafopy.interactive.funnel_chart(tf, width=WIDTH)
    .properties(title=f'Trim and fill: {int(numpy.sum(tf.fill))} studies added on the left')
    .configure_title(color=GRAPHITE, anchor='start')
)
funnel_h

hack = datasets.load('dat.hackshaw1998') fit_h = rma(yi='yi', vi='vi', data=hack) tf = trimfill(fit_h) funnel_h = ( metafopy.interactive.funnel_chart(tf, width=WIDTH) .properties(title=f'Trim and fill: {int(numpy.sum(tf.fill))} studies added on the left') .configure_title(color=GRAPHITE, anchor='start') ) funnel_h

Out[12]:

The usual tests and corrections side by side. The regression test (Egger et al., 1997, doi:10.1136/bmj.315.7109.629), the fail-safe N (Rosenthal, 1979) and a half-normal selection model (Vevea and Hedges, 1995; Preston et al., 2004) each look at the asymmetry in a different way.

In [13]:

summary_rows = pandas.DataFrame(
    [
        {'analysis': 'Random-effects model', 'estimate': float(fit_h.beta[0]), 'ci.lb': float(fit_h.ci_lb[0]), 'ci.ub': float(fit_h.ci_ub[0])},
        {'analysis': 'Trim and fill', 'estimate': float(tf.beta[0]), 'ci.lb': float(tf.ci_lb[0]), 'ci.ub': float(tf.ci_ub[0])},
        {
            'analysis': 'Selection model (half-normal)',
            **dict(
                zip(
                    ['estimate', 'ci.lb', 'ci.ub'],
                    [float(s) for s in (lambda sm: (sm.beta[0], sm.ci_lb[0], sm.ci_ub[0]))(selmodel(fit_h, type='halfnorm'))],
                    strict=True,
                )
            ),
        },
    ]
)
summary_rows[['estimate', 'ci.lb', 'ci.ub']] = numpy.exp(summary_rows[['estimate', 'ci.lb', 'ci.ub']])
summary_rows.columns = ['analysis', 'odds ratio', 'ci.lb', 'ci.ub']
summary_rows.round(3)

summary_rows = pandas.DataFrame( [ {'analysis': 'Random-effects model', 'estimate': float(fit_h.beta[0]), 'ci.lb': float(fit_h.ci_lb[0]), 'ci.ub': float(fit_h.ci_ub[0])}, {'analysis': 'Trim and fill', 'estimate': float(tf.beta[0]), 'ci.lb': float(tf.ci_lb[0]), 'ci.ub': float(tf.ci_ub[0])}, { 'analysis': 'Selection model (half-normal)', **dict( zip( ['estimate', 'ci.lb', 'ci.ub'], [float(s) for s in (lambda sm: (sm.beta[0], sm.ci_lb[0], sm.ci_ub[0]))(selmodel(fit_h, type='halfnorm'))], strict=True, ) ), }, ] ) summary_rows[['estimate', 'ci.lb', 'ci.ub']] = numpy.exp(summary_rows[['estimate', 'ci.lb', 'ci.ub']]) summary_rows.columns = ['analysis', 'odds ratio', 'ci.lb', 'ci.ub'] summary_rows.round(3)

Out[13]:

	analysis	odds ratio	ci.lb	ci.ub
0	Random-effects model	1.245	1.130	1.371
1	Trim and fill	1.191	1.083	1.309
2	Selection model (half-normal)	1.106	0.879	1.393

In [14]:

regtest(fit_h)

regtest(fit_h)

Out[14]:

Regression Test for Funnel Plot Asymmetry

Model:     mixed-effects meta-regression model
Predictor: standard error

Test for Funnel Plot Asymmetry: z = 2.1055, p = 0.0352
Limit Estimate (as sei -> 0):   b = 0.0216 (CI: -0.1774, 0.2207)

In [15]:

fsn(fit_h)

fsn(fit_h)

Out[15]:

Fail-safe N Calculation Using the General Approach

Average Effect Size:         0.2189 (with file drawer: 0.0473)
Amount of Heterogeneity:     0.0224 (with file drawer: 0.0312)
Observed Significance Level: <.0001 (with file drawer: 0.0504)
Target Significance Level:   0.05

Fail-safe N: 112

A generated write-up¶

metafor's reporter() writes a report of a fit; metafopy does not port it, because a notebook does the job better: the text below is built from the fits above, so it cannot disagree with them.

In [16]:

p = predict(fit, transf=numpy.exp)
q = mr.QM
qe_p = '< .0001' if fit.QEp < 1e-4 else f'= {fit.QEp:.4f}'
report = f"""
### Results

We pooled {fit.k} randomised trials of BCG vaccination with a random-effects model ({fit.method} estimate of τ²). The pooled risk ratio
of tuberculosis was {float(p.pred[0]):.2f} (95% CI {float(p.ci_lb[0]):.2f} to {float(p.ci_ub[0]):.2f}; 95% prediction interval
{float(p.pi_lb[0]):.2f} to {float(p.pi_ub[0]):.2f}). The trials were heterogeneous (Q({fit.k - 1}) = {fit.QE:.1f}, p {qe_p};
τ² = {fit.tau2:.3f}; I² = {fit.I2:.1f}%).

Absolute latitude explained a large part of the heterogeneity: each additional degree of latitude changed the log risk ratio by
{float(mr.beta[1]):.3f} (95% CI {float(mr.ci_lb[1]):.3f} to {float(mr.ci_ub[1]):.3f}; QM(1) = {q:.1f}), with R² = {mr.R2:.1f}%.

Across the {len(pooled['method'].unique())} τ² estimators the pooled risk ratio ranged from
{pooled.query("drop == '(none)'")['est'].min():.3f} to {pooled.query("drop == '(none)'")['est'].max():.3f}.
"""
IPython.display.Markdown(report)

p = predict(fit, transf=numpy.exp) q = mr.QM qe_p = '< .0001' if fit.QEp < 1e-4 else f'= {fit.QEp:.4f}' report = f""" ### Results We pooled {fit.k} randomised trials of BCG vaccination with a random-effects model ({fit.method} estimate of τ²). The pooled risk ratio of tuberculosis was {float(p.pred[0]):.2f} (95% CI {float(p.ci_lb[0]):.2f} to {float(p.ci_ub[0]):.2f}; 95% prediction interval {float(p.pi_lb[0]):.2f} to {float(p.pi_ub[0]):.2f}). The trials were heterogeneous (Q({fit.k - 1}) = {fit.QE:.1f}, p {qe_p}; τ² = {fit.tau2:.3f}; I² = {fit.I2:.1f}%). Absolute latitude explained a large part of the heterogeneity: each additional degree of latitude changed the log risk ratio by {float(mr.beta[1]):.3f} (95% CI {float(mr.ci_lb[1]):.3f} to {float(mr.ci_ub[1]):.3f}; QM(1) = {q:.1f}), with R² = {mr.R2:.1f}%. Across the {len(pooled['method'].unique())} τ² estimators the pooled risk ratio ranged from {pooled.query("drop == '(none)'")['est'].min():.3f} to {pooled.query("drop == '(none)'")['est'].max():.3f}. """ IPython.display.Markdown(report)

Out[16]:

Results¶

We pooled 13 randomised trials of BCG vaccination with a random-effects model (REML estimate of τ²). The pooled risk ratio of tuberculosis was 0.49 (95% CI 0.34 to 0.70; 95% prediction interval 0.15 to 1.55). The trials were heterogeneous (Q(12) = 152.2, p < .0001; τ² = 0.313; I² = 92.2%).

Absolute latitude explained a large part of the heterogeneity: each additional degree of latitude changed the log risk ratio by -0.029 (95% CI -0.043 to -0.015; QM(1) = 16.4), with R² = 75.6%.

Across the 7 τ² estimators the pooled risk ratio ranged from 0.488 to 0.491.

Exported figures¶

The same fits drive static figures for a paper: the matplotlib forest plot of metafopy, and the Altair charts saved as SVG or PNG with vl-convert (no browser needed). The files go to a temporary directory; point out anywhere you like.

In [17]:

import matplotlib
import matplotlib.pyplot

from metafopy import forest

out = pathlib.Path(tempfile.mkdtemp(prefix='metafopy-showcase-'))
fit_labelled = rma(yi='yi', vi='vi', data=es, slab=es['study'])
matplotlib.pyplot.figure(figsize=(8, 6))
forest(fit_labelled, atransf=numpy.exp, header=True)
matplotlib.pyplot.savefig(out / 'forest.png', dpi=150, bbox_inches='tight')
matplotlib.pyplot.close()

dashboard.save(out / 'dashboard.svg')
dashboard.save(out / 'dashboard.png', scale_factor=2)
funnel_h.save(out / 'funnel-hackshaw.svg')
print(*sorted(f'{f.name} ({f.stat().st_size / 1024:.0f} KB)' for f in out.iterdir()), sep='\n')
IPython.display.Image(out / 'forest.png', width=520)

import matplotlib import matplotlib.pyplot from metafopy import forest out = pathlib.Path(tempfile.mkdtemp(prefix='metafopy-showcase-')) fit_labelled = rma(yi='yi', vi='vi', data=es, slab=es['study']) matplotlib.pyplot.figure(figsize=(8, 6)) forest(fit_labelled, atransf=numpy.exp, header=True) matplotlib.pyplot.savefig(out / 'forest.png', dpi=150, bbox_inches='tight') matplotlib.pyplot.close() dashboard.save(out / 'dashboard.svg') dashboard.save(out / 'dashboard.png', scale_factor=2) funnel_h.save(out / 'funnel-hackshaw.svg') print(*sorted(f'{f.name} ({f.stat().st_size / 1024:.0f} KB)' for f in out.iterdir()), sep='\n') IPython.display.Image(out / 'forest.png', width=520)

dashboard.png (578 KB)
dashboard.svg (101 KB)
forest.png (108 KB)
funnel-hackshaw.svg (36 KB)

Out[17]:

Where next¶

• User guide and function reference for everything not shown here.
• Coming from R if you know metafor already.
• Run this notebook yourself: pip install "metafopy[interactive] @ git+https://github.com/bnomis/metafopy" jupyterlab, or use the Binder badge at the top.