.. _confirmative:

Population inference: Confirmation study
----------------------------------------




The analysis in :ref:`subject-inference` suggested that a word
generating task results in a relevant and statistical significant
increase in BOLD effect in the occipital lobe:


.. image:: figures/population-inference-confirmative_figure2_1.png
   :width: 15 cm



In this section, we aim to confirm this finding in a population of
right-handed, healthy subjects. For right-handed, healthy subjects (with
a `Waterloo Handedness Score`_ above +20) from the population have been
drawn, their FMRI have been analysed, and the results collected in the
file ``tutorial.sample``:

.. code::

    fmrisample -v tutorial.sample \
        --covariates-query 'waterloo > 20'

.. _`Waterloo Handedness Score`: https://www.ucl.ac.uk/medical-education/resources/Waterloo/WatFoot_HandQuest36items-Elias1998.pdf

Load data and fit the population model
......................................


.. code:: python

    from fmristats import *
    from fmristats.plot import *
    
    population_sample = load(path_to_sample)
    
    print("""
        Number of subjects in total:     {:d}
        Number of right-handed subjects: {:d}
        Number of  left-handed subjects: {:d}""".format(
            len(population_sample.covariates),
            (population_sample.covariates.waterloo >= 20).sum(),
            (population_sample.covariates.waterloo <= 20).sum()))


.. code::

    
        Number of subjects in total:     64
        Number of right-handed subjects: 64
        Number of  left-handed subjects: 0
    



We are aiming to confirm a finding that has been suggested by an
analysis of a single subject. Hence, remove this *initial subject* from
the sample.


.. code:: python

    
    sample = population_sample.filter(population_sample.covariates.id != 2)




First, we shall have a look at the sample size per voxel in the
MNI152-template (the number of subjects which have valid entries at this
voxel):


.. code:: python

    
    sample_size = Image(
        sample.vb.reference,
        np.isfinite(sample.statistics).all(axis=-2).sum(axis=-1))
    
    picture(sample_size.mask(inplace=False),cmap=cm.viridis)


.. image:: figures/population-inference-confirmative_figure5_1.png
   :width: 15 cm



It appears that for some voxels close to the surface of the
MNI152-template the population diffeomorphisms :math:`ψ_j` did not fully
map the template brain onto the respective subject brain. This is to be
expected.


.. code:: python

    
    print('Sample size at the index of interest: {:.0f}'.format(
        sample_size.data[index]))


.. code::

    Sample size at the index of interest: 63
    



Let us now fit the population model (which is a random effects meta
regression model in fmristats) to the sample of BOLD effects and their
respective estimated standard errors:


.. code:: python

    
    popmod = PopulationModel(sample)
    result = popmod.fit()


.. code::

      … not all points in the mask are identifiable.
      … points with missing data along subject dimension.
      … number of points to estimate: 902629
      … perform a meta analysis
    



Extract the Knapp-Hartung adjusted t-statistics field
.....................................................

Extract the Knapp-Hartung adjusted t-statistics field that tests for
non-zero, task related BOLD effects:


.. code:: python

    tstatistics = result.get_tstatistic()
    
    picture(tstatistics)


.. image:: figures/population-inference-confirmative_figure8_1.png
   :width: 15 cm



Now, plot the respective slice through the t-statistics field that has
shown the above BOLD effect in the initial subject:


.. code:: python

    
    pt.plot([index[0]], [index[1]], 'ko')
    pt.plot([index[0]], [index[1]], 'w+')
    picture(tstatistics,3,1,1,[index[-1]],
        interpolation='bilinear')


.. image:: figures/population-inference-confirmative_figure9_1.png
   :width: 15 cm



The t-statistic at these coordinates is:


.. code:: python

    print(tstatistics.data[index])


.. code::

    2.705125720068511
    






Create a forest plot of BOLD statistics
.......................................

Extract the estimated BOLD effects and the respective standard errors of
the estimated BOLD effects from the sample:


.. code:: python

    
    df = population_sample.at_index(index)
    df.valid.all()
    df.sort_values(by='waterloo', inplace=True)
    
    dm = result.at_index(index)
    dm.set_index('parameter', inplace=True)




Define the critical values for the plots:


.. code:: python

    
    from scipy.stats.distributions import t
    from scipy.stats.distributions import norm
    
    crt_subject = norm.ppf(q=.95)
    crt_population = t.ppf(q=.95, df = dm.loc['Intercept', 'df'])
    
    print(crt_subject, crt_population)


.. code::

    1.6448536269514722 1.669804162296528
    



Set some general options for the visualisation:


.. code:: python

    
    import matplotlib
    matplotlib.rc('xtick', labelsize=8)
    matplotlib.rc('ytick', labelsize=8)
    
    import seaborn as sb
    sb.set_style('whitegrid')
    palette = sb.palettes.SEABORN_PALETTES['deep']
    
    figw = 5.842
    figh = 8.442




Create a forest plot in ascending order of handiness (left handed on the
bottom, right handed on the top):


.. code:: python

    
    x = dm.loc['Intercept', 'point']
    s = dm.loc['Intercept', 'stderr']
    
    df['yvec'] = range(len(df.task))
    df['waterloo'] = df.waterloo.astype(int)
    
    fig = pt.figure(figsize=(figw,figh))
    ax = pt.subplot(111)
    ax.axvline(0,c='k',lw=.9, ls=':')
    ax.errorbar(df.task[df.id!=2], df.yvec[df.id!=2],
        xerr=crt_subject*df.stderr[df.id!=2], fmt='o',
        label='Sample', c=palette[0])
    ax.errorbar(df.task[df.id==2], df.yvec[df.id==2],
        xerr=crt_subject*df.stderr[df.id==2], fmt='o',
        label='Initial Subject', c=palette[1])
    ax.errorbar(x, -1, xerr=crt_population*s, fmt='o',
        label='Population', c=palette[2])
    ax.axhline(-.5,c='k',lw=.9, ls='-')
    ax.set_xlabel(r'Task related BOLD effect at $x_0$ in ATI')
    ax.set_ylabel('Waterloo')
    ax.yaxis.set_ticks_position('none')
    pt.box(False)
    pt.yticks(np.hstack((df.yvec, -1)), list(df.waterloo) + ["Meta"] )
    pt.legend(loc='lower center', bbox_to_anchor=(.5, -0.2), ncol=3)


.. image:: figures/population-inference-confirmative_figure15_1.png
   :width: 15 cm






Quite interestingly, there was even one subject in the sample that
showed an even larger task related BOLD effect than the initial subject
at theses coordinates.

Create a funnel plot of BOLD statistics
.......................................

Create a funnel plot with the estimated BOLD effect on the x-axis and
the reciprocal of the squared standard error on the y-axis:





.. image:: figures/population-inference-confirmative_figure18_1.png
   :width: 15 cm



Formal testing
..............

Any prior knowledge had been removed prior to testing the hypothesis of
whether there exists a task related BOLD effect at :math:`x` in the
sample. Hence, there is no need to adjust for multiple testing in the
inference for the population, we there was, well, no multiple testing.

Point estimate, confidence interval, and p-value of BOLD effect of the
word generating task at the point :math:`x` in the population:


.. code:: python

    x = dm.loc['Intercept', 'point']
    s = dm.loc['Intercept', 'stderr']
    
    print("""
    Point estimate of the BOLD effect: {:.2f}
    Lower bound of a 95% confidence interval: {:.2f}
    Upper bound of a 95% confidence interval: {:.2f}
    p-value for ≠0:                           {:.4f}""".format(
        x, x-crt_population*s, x+crt_population*s,
        dm.loc['Intercept', 'pvalue']))


.. code::

    
    Point estimate of the BOLD effect: 5.97
    Lower bound of a 95% confidence interval: 2.29
    Upper bound of a 95% confidence interval: 9.66
    p-value for ≠0:                           0.0044
    



We may reject the null hypothesis that there exists no BOLD effect at
:math:`x` at :math:`1\%` level of significance.