# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:python3]
#     language: python
#     name: conda-env-python3-py
# ---

# # Run meta-analyses on 21 pain studies
# Collection of NIDM-Results packs downloaded from Neurovault collection [1425](https://www.neurovault.org/collections/1425/), uploaded by Dr. <NAME>.

# +
# %matplotlib inline
import json
import numpy as np
from glob import glob
from os.path import basename, join, dirname, isfile

import pandas as pd
import nibabel as nib
from scipy.stats import t
from nilearn.masking import apply_mask
from nilearn.plotting import plot_stat_map

import nimare
from nimare.meta.ibma import (stouffers, fishers, weighted_stouffers,
                              rfx_glm, ffx_glm)

# +
dset_file = '/Users/tsalo/Documents/tsalo/NiMARE/nimare/tests/data/nidm_pain_dset.json'
with open(dset_file, 'r') as fo:
    dset_dict = json.load(fo)
db = nimare.dataset.Database(dset_file)
dset = db.get_dataset()

mask_img = dset.mask

logp_thresh = -np.log(.05)


# -

# ## Define functions for querying set of NIDM-Results packs

# +
def _get_file(cdict, t):
    """Return the file associated with a given data type within a
    folder if it exists. Otherwise, returns an empty list.
    """
    temp = ''
    if t == 'con':
        temp = cdict['images'].get('con')
    elif t == 'se':
        temp = cdict['images'].get('se')
    elif t == 't':
        temp = cdict['images'].get('t')
    elif t == 'z':
        temp = cdict['images'].get('z')
    elif t == 'n':
        temp = cdict.get('sample_sizes', [])
        if temp:
            temp = np.mean(temp)
    else:
        raise Exception('Input type "{0}" not recognized.'.format(t))

    return temp


def get_files(ddict, types):
    """Returns a list of files associated with a given data type 
    from a set of subfolders within a directory. Allows for
    multiple data types and only returns a set of files from folders
    with all of the requested types.
    """
    all_files = []
    for study in ddict.keys():
        files = []
        cdict = ddict[study]['contrasts']['1']
        for t in types:
            temp = _get_file(cdict, t)
            if temp:
                files.append(temp)
            
        if len(files) == len(types):
            all_files.append(files)
    all_files = list(map(list, zip(*all_files)))
    return all_files


# -

# ## Get z-maps

files = get_files(dset_dict, ['z'])
z_imgs = [nib.load(f) for f in files[0]]
z_data = apply_mask(z_imgs, mask_img)
print('{0} studies found.'.format(z_data.shape[0]))

# ## Fisher's

# %%time
result = fishers(z_data, mask_img)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ## Stouffer's

# ### Fixed-effects inference

# %%time
result = stouffers(z_data, mask_img, inference='ffx',
                   null='theoretical', n_iters=None)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ### Random-effects inference with theoretical null distribution

# %%time
result = stouffers(z_data, mask_img, inference='rfx',
                   null='theoretical', n_iters=None)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ### Random-effects inference with empirical null distribution

# %%time
result = stouffers(z_data, mask_img, inference='rfx',
                   null='empirical', n_iters=1000)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ## Get z-maps + sample sizes

z_files, ns = get_files(dset_dict, ['z', 'n'])
z_imgs = [nib.load(f) for f in z_files]
z_data = apply_mask(z_imgs, mask_img)
sample_sizes = np.array(ns)
print('{0} studies found.'.format(z_data.shape[0]))

# ## Weighted Stouffer's

# %%time
result = weighted_stouffers(z_data, sample_sizes, mask_img)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ## Get contrast maps + contrast standard error maps + sample sizes

con_files, se_files, ns = get_files(dset_dict, ['con', 'se', 'n'])
con_imgs = [nib.load(f) for f in con_files]
se_imgs = [nib.load(f) for f in se_files]
con_data = apply_mask(con_imgs, mask_img)
se_data = apply_mask(se_imgs, mask_img)
sample_sizes = np.array(ns)
print('{0} studies found.'.format(con_data.shape[0]))

# ## FFX GLM

# %%time
result = ffx_glm(con_data, se_data, sample_sizes, mask_img, equal_var=True)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ## Get contrast maps

con_files = get_files(dset_dict, ['con'])
con_files = con_files[0]
con_imgs = [nib.load(f) for f in con_files]
con_data = apply_mask(con_imgs, mask_img)
print('{0} studies found.'.format(con_data.shape[0]))

# ## RFX GLM

# ### Theoretical null distribution

# %%time
result = rfx_glm(con_data, mask_img, null='theoretical', n_iters=None)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ### Empirical null distribution

# %%time
result = rfx_glm(con_data, mask_img, null='empirical', n_iters=1000)
plot_stat_map(result.images['log_p'], threshold=logp_thresh,
              cut_coords=[0, 0, -8], draw_cross=False,
              cmap='RdBu_r')

# ## MKDA

# ### Density analysis

# %%time
mkda = nimare.meta.cbma.MKDADensity(dset, ids=dset.ids, kernel__r=10)
mkda.fit(n_iters=10)
plot_stat_map(mkda.results.images['vfwe'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

# ### Chi2 with FDR correction

# %%time
mkda2 = nimare.meta.cbma.MKDAChi2(dset, ids=dset.ids, ids2=dset.ids, kernel__r=10)
mkda2.fit(corr='FDR')
plot_stat_map(mkda2.results.images['pAgF_z_FDR_05'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

# ### Chi2 with FWE correction

# %%time
mkda3 = nimare.meta.cbma.MKDAChi2(dset, ids=dset.ids, ids2=dset.ids, kernel__r=10)
mkda3.fit(corr='FWE', n_iters=10)
plot_stat_map(mkda3.results.images['pAgF_z_FWE_05'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

# ## KDA

# %%time
kda = nimare.meta.cbma.KDA(dset, ids=dset.ids, kernel__r=10)
kda.fit(n_iters=10)
plot_stat_map(kda.results.images['vfwe'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

# ## ALE

# %%time
ale = nimare.meta.cbma.ALE(dset, ids=dset.ids)
ale.fit(n_iters=10)
plot_stat_map(ale.results.images['vfwe'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

# ## SCALE

# %%time
ijk = np.vstack(np.where(dset.mask.get_data())).T
scale = nimare.meta.cbma.SCALE(dset, ids=dset.ids, ijk=ijk)
scale.fit(n_iters=10)
plot_stat_map(scale.results.images['vthresh'], cut_coords=[0, 0, -8],
              draw_cross=False, cmap='RdBu_r')

StarcoderdataJupyterScriptsDedupFiltered