Quantitative comparisons¶

Can we express quantitatively what we observed in the previous plots? This is where the meta-analysis tools come in: we used the R package metafor to fit a mixed-effect (ME) model to the data reported for each measure. In this way, we can estimate an overall interval of R² values based on the effect sizes and the sample sizes. We can also estimate the interval of R² that we can expect in future studies (this is called prediction interval).

Figure 5

A compact way to represent these results is given by forest plots: for each study, we represent the effect size and the related sample size using a square and a horizontal error bar; then for each measure, we represent the results from the ME model using a diamond and an additional error bar; finally to represent the prediction interval we use two hourglasses and a dotted line.

import numpy as np
import pandas as pd

import plotly.graph_objects as go
from IPython.core.display import display, HTML
from plotly.offline import plot
import plotly.express as px
import plotly.colors
from plotly.subplots import make_subplots

from rpy2.robjects.packages import importr
import rpy2.robjects
import subprocess
subprocess.call('curl https://raw.githubusercontent.com/Notebook-Factory/brand/main/insertLogo.py --output /tmp/insertLogo.py', shell=True)
%run /tmp/insertLogo.py

Figure 5¶

config={'showLink': False, 'displayModeBar': False}

filtered_df = pd.read_pickle("filtered_df.pkl")

filtered_df['Variance'] = (4*filtered_df['R^2'])*((1-filtered_df['R^2'])**2)/filtered_df['Sample points']

measure_type = {'Diffusion':['RD', 'AD', 'FA', 'MD',
                'AWF', 'RK', 'RDe', 'MK'],
                'Magnetization transfer':['MTR',
                'ihMTR', 'MTR-UTE', 'MPF', 'MVF-MT',
                'R1f', 'T2m', 'T2f', 'k_mf','k_fm'],
                'T1 relaxometry':['T1'], 'T2 relaxometry':['T2', 'MWF', 'MVF-T2'],
                'Other':['QSM', 'R2*', 'rSPF', 'MTV',
                'T1p', 'T2p', 'RAFF', 'PD', 'T1sat']}

color_dict = {m:plotly.colors.qualitative.Bold[n]
              for n,m in enumerate(measure_type.keys())}

metafor = importr('metafor')
stats = importr('stats')

metastudy = {}
for m in filtered_df.Measure.unique():
    nstudies=len(filtered_df.Measure[filtered_df.Measure==m])
    if nstudies > 2:
        df_m = filtered_df[filtered_df.Measure==m]
        df_m = df_m.sort_values(by=['Year'])
        
        r2 = rpy2.robjects.FloatVector(df_m['R^2'])
        var = rpy2.robjects.FloatVector(df_m['Variance'])
        fit = metafor.rma(r2, var, method="REML", test="knha")
        res = stats.predict(fit)
        
        results = dict(zip(res.names,list(res)))
        
        metastudy[m] = dict(pred=results['pred'][0], cilb=results['pred'][0]-results['ci.lb'][0],
                            ciub=results['ci.ub'][0]-results['pred'][0],
                            crub=results['cr.ub'][0], 
                            crlb=results['cr.lb'][0])
measure_type_reverse={m:t for t,mlist in measure_type.items() for m in mlist}
fig5 = make_subplots(rows=3, cols=3, start_cell="top-left", vertical_spacing=0.06,
                      horizontal_spacing=0.2, x_title='R<sup>2</sup>',
                      subplot_titles=sorted(metastudy.keys(), key=measure_type_reverse.get))

row=1
col=1
for m in sorted(metastudy.keys(), key=measure_type_reverse.get):
    fig5.add_trace(go.Scatter(
        x=[round(metastudy[m]['crlb'],2) if round(metastudy[m]['crlb'],2)>0 else 0,
           round(metastudy[m]['crub'],2) if round(metastudy[m]['crub'],2)<1 else 1],
        y=['Mixed model','Mixed model'],
        line=dict(color='black', width=2, dash='dot'),
        hovertemplate = 'Prediction boundary: %{x}<extra></extra>',
        marker_symbol = 'hourglass-open', marker_size = 8
    ), row=row, col=col)
    
    fig5.add_trace(go.Scatter(
        x=[round(metastudy[m]['pred'],2)],
        y=['Mixed model'],
        mode='markers',
        marker = dict(color = 'black'),
        marker_symbol = 'diamond-wide',
        marker_size = 10,
        hovertemplate = 'R<sup>2</sup> estimate: %{x}<extra></extra>',
        error_x=dict(
            type='data',
            arrayminus=[round(metastudy[m]['cilb'],2) if round(metastudy[m]['cilb'],2)>0 else 0],
            array=[round(metastudy[m]['ciub'],2) if round(metastudy[m]['ciub'],2)<1 else 1])
    ), row=row, col=col)
    
    df_m = filtered_df[filtered_df.Measure==m]
    df_m = df_m.sort_values(by=['Year'], ascending=False)
    fig5.add_trace(go.Scatter(
        x=df_m['R^2'],
        y=df_m['Study'],
        text=df_m['Sample points'],
        customdata=df_m['Histology/microscopy measure'],
        mode='markers',
        marker = dict(color = color_dict[measure_type_reverse[m]]),
        marker_symbol = 'square',
        marker_size = np.log(50/df_m['Variance']),
        hovertemplate = '%{y}<br>R<sup>2</sup>: %{x}<br>Number of samples: %{text}' +
            '<br>Reference: %{customdata}<extra></extra>',
        error_x=dict(
            type='data',
            array=2*np.sqrt(df_m['Variance']))
    ), row=row, col=col)

    if col == 3:
        col = 1
        row += 1
    else:
        col += 1

fig5.update_xaxes(range=[0, 1])

fig5.update_yaxes(tickfont=dict(size=8))

fig5.update_layout(showlegend=False,
    title=dict(text='Figure 5: Forest plots and mixed modelling results'),
    margin=dict(l=0),
    width=700,
    height=800)

for ii in range(9):
    fig5.layout.annotations[ii]["font"] = {'size': 10,'color':'black'}

plot(insertLogo(fig5,0.03,0.03,1,-0.07,-0.125,0.033), filename = 'fig5.html',config = config)
display(HTML('fig5.html'))

The quest for measuring myelin with MRI - An interactive meta-analysis

Quantitative comparisons¶

Figure 5¶