{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nTutorial 3: Null models for gradient significance\n==================================================\nIn this tutorial we assess the significance of correlations between the first\ncanonical gradient and data from other modalities (curvature, cortical\nthickness and T1w/T2w image intensity). A normal test of the significance of\nthe correlation cannot be used, because the spatial auto-correlation in MRI\ndata may bias the test statistic. In this tutorial we will show two approaches\nfor null hypothesis testing: spin permutations and Moran spectral\nrandomization.\n\n

Note

When using either approach to compare gradients to non-gradient markers,\n we recommend randomizing the non-gradient markers as these randomizations\n need not maintain the statistical independence between gradients.

\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spin Permutations\n------------------------------\n\nHere, we use the spin permutations approach previously proposed in\n`(Alexander-Bloch et al., 2018)\n`_,\nwhich preserves the auto-correlation of the permuted feature(s) by rotating\nthe feature data on the spherical domain.\nWe will start by loading the conte69 surfaces for left and right hemispheres,\ntheir corresponding spheres, midline mask, and t1w/t2w intensity as well as\ncortical thickness data, and a template functional gradient.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom brainspace.datasets import load_gradient, load_marker, load_conte69\n\n# load the conte69 hemisphere surfaces and spheres\nsurf_lh, surf_rh = load_conte69()\nsphere_lh, sphere_rh = load_conte69(as_sphere=True)\n\n# Load the data\nt1wt2w_lh, t1wt2w_rh = load_marker('t1wt2w')\nt1wt2w = np.concatenate([t1wt2w_lh, t1wt2w_rh])\n\nthickness_lh, thickness_rh = load_marker('thickness')\nthickness = np.concatenate([thickness_lh, thickness_rh])\n\n# Template functional gradient\nembedding = load_gradient('fc', idx=0, join=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let\u2019s first generate some null data using spintest.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.null_models import SpinPermutations\nfrom brainspace.plotting import plot_hemispheres\n\n# Let's create some rotations\nn_rand = 1000\n\nsp = SpinPermutations(n_rep=n_rand, random_state=0)\nsp.fit(sphere_lh, points_rh=sphere_rh)\n\nt1wt2w_rotated = np.hstack(sp.randomize(t1wt2w_lh, t1wt2w_rh))\nthickness_rotated = np.hstack(sp.randomize(thickness_lh, thickness_rh))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an illustration of the rotation, let\u2019s plot the original t1w/t2w data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plot original data\nplot_hemispheres(surf_lh, surf_rh, array_name=t1wt2w, size=(1200, 200), cmap='viridis',\n nan_color=(0.5, 0.5, 0.5, 1), color_bar=True, zoom=1.65)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "as well as a few rotated versions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = 2\n# Plot some rotations\nplot_hemispheres(surf_lh, surf_rh, array_name=t1wt2w_rotated[:3], size=(1200, 600),\n cmap='viridis', nan_color=(0.5, 0.5, 0.5, 1), color_bar=True,\n zoom=1.55, label_text=['Rot0', 'Rot1', 'Rot2'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Warning

With spin permutations, midline vertices (i.e,, NaNs) from both the\n original and rotated data are discarded. Depending on the overlap of\n midlines in the, statistical comparisons between them may compare\n different numbers of features. This can bias your test statistics.\n Therefore, if a large portion of the sphere is not used, we recommend\n using Moran spectral randomization instead.

\n\nNow we simply compute the correlations between the first gradient and the\noriginal data, as well as all rotated data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\nfrom scipy.stats import spearmanr\n\nfig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w, 'thickness': thickness}\nrotated = {'t1wt2w': t1wt2w_rotated, 'thickness': thickness_rotated}\n\nr_spin = np.empty(n_rand)\nmask = ~np.isnan(thickness)\nfor k, (fn, feat) in enumerate(feats.items()):\n r_obs, pv_obs = spearmanr(feat[mask], embedding[mask])\n\n # Compute perm pval\n for i, perm in enumerate(rotated[fn]):\n mask_rot = mask & ~np.isnan(perm) # Remove midline\n r_spin[i] = spearmanr(perm[mask_rot], embedding[mask_rot])[0]\n pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))\n\n # Plot null dist\n axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))\n axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n axs[k].set_xlabel('Correlation with {}'.format(fn))\n if k == 0:\n axs[k].set_ylabel('Density')\n\n print('{}:\\n Obs : {:.5e}\\n Spin: {:.5e}\\n'.\n format(fn.capitalize(), pv_obs, pv_spin))\n\nfig.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is interesting to see that both p-values increase when taking into\nconsideration the auto-correlation present in the surfaces. Also, we can see\nthat the correlation with thickness is no longer statistically significant\nafter spin permutations.\n\n\n\nMoran Spectral Randomization\n------------------------------\n\nMoran Spectral Randomization (MSR) computes Moran's I, a metric for spatial\nauto-correlation and generates normally distributed data with similar\nauto-correlation. MSR relies on a weight matrix denoting the spatial\nproximity of features to one another. Within neuroimaging, one\nstraightforward example of this is inverse geodesic distance i.e. distance\nalong the cortical surface.\n\nIn this example we will show how to use MSR to assess statistical\nsignificance between cortical markers (here curvature and cortical t1wt2w\nintensity) and the first functional connectivity gradient. We will start by\nloading the left temporal lobe mask, t1w/t2w intensity as well as cortical\nthickness data, and a template functional gradient\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.datasets import load_mask\n\nn_pts_lh = surf_lh.n_points\nmask_tl, _ = load_mask(name='temporal')\n\n# Keep only the temporal lobe.\nembedding_tl = embedding[:n_pts_lh][mask_tl]\nt1wt2w_tl = t1wt2w_lh[mask_tl]\ncurv_tl = load_marker('curvature')[0][mask_tl]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now compute the Moran eigenvectors. This can be done either by\nproviding a weight matrix of spatial proximity between each vertex, or by\nproviding a cortical surface. Here we\u2019ll use a cortical surface.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from brainspace.null_models import MoranRandomization\nfrom brainspace.mesh import mesh_elements as me\n\n# compute spatial weight matrix\nw = me.get_ring_distance(surf_lh, n_ring=1, mask=mask_tl)\nw.data **= -1\n\n\nmsr = MoranRandomization(n_rep=n_rand, procedure='singleton', tol=1e-6,\n random_state=0)\nmsr.fit(w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the Moran eigenvectors we can now compute the randomized data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "curv_rand = msr.randomize(curv_tl)\nt1wt2w_rand = msr.randomize(t1wt2w_tl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the randomized data, we can compute correlations between\nthe gradient and the real/randomised data and generate the non-parametric\np-values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}\nrand = {'t1wt2w': t1wt2w_rand, 'curvature': curv_rand}\n\nfor k, (fn, data) in enumerate(rand.items()):\n r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')\n\n # Compute perm pval\n r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])\n pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))\n\n # Plot null dist\n axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))\n axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n axs[k].set_xlabel('Correlation with {}'.format(fn))\n if k == 0:\n axs[k].set_ylabel('Density')\n\n print('{}:\\n Obs : {:.5e}\\n Moran: {:.5e}\\n'.\n format(fn.capitalize(), pv_obs, pv_rand))\n\nfig.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are some scenarios where MSR results do not follow a normal\ndistribution. It is relatively simple to check whether this occurs in our\ndata by visualizing the null distributions. Check this interesting paper\nfor more information `(Burt et al., 2020) `_.\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }