Quantitative comparisons
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 R2 values based on the effect sizes and the sample sizes. We can also estimate the interval of R2 that we can expect in future studies (this is called prediction interval).
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
import os
from repo2data.repo2data import Repo2Data
data_req_path = os.path.join("..","..", "binder", "data_requirement.json")
repo2data = Repo2Data(data_req_path)
data_path = repo2data.install()[0]
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.py1Figure 5¶
config={'showLink': False, 'displayModeBar': False}
filtered_df = pd.read_pickle(os.path.join(data_path,'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'))Loading...
The forest plot offers a detailed summary for each measure. What if we want to compare the R2 estimates across measures? To do that, we pooled together all the measures from all the studies and computed first a repeated measures meta-regression and then all the possible pairwise comparisons (Tukey’s test), correcting for multiple comparisons (Bonferroni correction).
2Figure 6¶
multcomp = importr('multcomp')
base = importr('base')
thres = filtered_df.Measure.value_counts() > 2
df_thres = filtered_df[filtered_df.Measure.isin(filtered_df.Measure.value_counts()[thres].index)]
r2 = rpy2.robjects.FloatVector(df_thres['R^2'])
var = rpy2.robjects.FloatVector(df_thres['Variance'])
measure_v = rpy2.robjects.StrVector(df_thres['Measure'])
measure_f = rpy2.robjects.Formula('~ -1 + measure')
env = measure_f.environment
env['measure'] = measure_v
study_v = rpy2.robjects.StrVector(df_thres['Study'])
study_f = rpy2.robjects.Formula('~ 1 | study')
env = study_f.environment
env['study'] = study_v
fit_mv = metafor.rma_mv(r2, var, method="REML", mods=measure_f, random=study_f)
glht = multcomp.glht(fit_mv, base.cbind(multcomp.contrMat(base.rep(1,9), type="Tukey")))
mtests = multcomp.summary_glht(glht, test=multcomp.adjusted("bonferroni"))
mtest_res = dict(zip(mtests.names, list(mtests)))
stat_res = dict(zip(mtest_res['test'].names, list(mtest_res['test'])))
pvals = stat_res['pvalues']
zvals = stat_res['tstat']
measure_list = df_thres['Measure'].unique()
measure_list.sort()
n = len(measure_list)
pvals_list = []
zvals_list = []
idx = 0
for i in range(n):
pvals_i = [0] * i
pvals_i.append(np.nan)
pvals_i.extend(pvals[idx:idx+n-i-1])
zvals_i = [0] * i
zvals_i.append(np.nan)
zvals_i.extend(zvals[idx:idx+n-i-1])
idx = idx + n - i - 1
pvals_list.append(pvals_i)
zvals_list.append(zvals_i)
pvals_list=np.array(pvals_list)
zvals_list=np.array(zvals_list)
pvals_list=pvals_list+pvals_list.T
zvals_list=zvals_list-zvals_list.T
fig6 = make_subplots(rows=1, cols=2, horizontal_spacing=0.14,
subplot_titles=['z-scores',
'p-values'
])
fig6.add_trace(go.Heatmap(
z=zvals_list,
x=measure_list,
y=measure_list,
customdata=pvals_list,
hovertemplate = '%{x}-%{y}<br>z-score: %{z:.2f}<br>p-value: %{customdata:.5f}<extra></extra>',
hoverongaps = False,
colorscale='RdBu',
colorbar=dict(title='z-score', x=0.42,thickness=20),
showscale=True
), col=1, row=1)
fig6.add_trace(go.Heatmap(
z=pvals_list,
x=measure_list,
y=measure_list,
customdata=zvals_list,
hovertemplate = '%{x}-%{y}<br>z-score: %{customdata:.2f}<br>p-value: %{z:.5f}<extra></extra>',
hoverongaps = False,
colorscale='Purples',
colorbar=dict(title='p-value',thickness=20),
zmin=0,
zmax=0.05,
showscale=True,
reversescale=True
), col=2, row=1)
fig6.update_layout(
title='Figure 6: Statistical pairwise comparisons between R<sup>2</sup> estimates',
width=740,
height = 420,
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(0,0,0,0)',
margin=dict(l=0,r=0)
)
plot(insertLogo(fig6,0.07,0.07,1.11,-0.15,-0.1,0.075), filename = 'fig6.html',config = config)
display(HTML('fig6.html'))Loading...