notebooks/nb05_Group-Level-Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Author:** [Dace Apšvalka](https://www.mrc-cbu.cam.ac.uk/people/dace.apsvalka/) with edits by [Rik Henson](https://www.mrc-cbu.cam.ac.uk/people/rik.henson/)\n",
    "- **Date:** February 2025  \n",
    "- **conda environment**: [mri conda environment](https://github.com/dcdace/fMRI_training/blob/main/mri_environment.yml)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fMRI Data Analysis: Group-Level Analysis\n",
    "\n",
    "Once you have beta (or contrast) maps for conditions (or contrasts) from all subjects, you can perform group-level statistics. Importantly, all subject first-level results need to be in common space, e.g., MNI, to perform voxel-wise group analyses. Group-level analysis allows you to make inferences about the population, rather than individual subjects, by assessing common activations across participants. Common statistical methods for group-level analysis include one-sample or paired t-tests, as well as more complex ANOVA models, depending on your study design.\n",
    "\n",
    "In this tutorial, we illustrate an ANOVA in which we include 9 contrasts, one for each of the 9 experimental conditions (averaged across runs), as well as 16 additional regressors to remove between-subject variance. This corresponds to a *repeated-measures ANOVA*, though assumes that the (pooled) error is spherical *(see [Rik's Stats tutorial](https://github.com/RikHenson/CognesticStats) section on ANOVA for the importance of this)*. \n",
    "\n",
    "Here is a recommended viewing to help better understand the principles of the group-level analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAACAwEBAQAAAAAAAAAAAAAAAwECBAUGB//EADoQAAIBAwEGBgEDAwIEBwAAAAABAgMREgQFEyExUZEVIkFSYdEUMkJxBoGhI2IkM7HwFkNyksHh8f/EABkBAQEBAQEBAAAAAAAAAAAAAAABAgMEBf/EACQRAQEBAQADAQACAgIDAAAAAAABEQIDITESQVEiMnGRBBNh/9oADAMBAAIRAxEAPwD5+AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA/8WfWJP4dTrHuBnA0/hVesO5P4FX3Q7sDKBq8Pq+6Hdk+HVvdDu/oDIBs8Ore6Hd/QeG1vdDu/oYMYG3wyt7qfd/QeGVvdT7v6LgxAbfDK/up939B4XX91Pu/oYMQG3wuv7qfd/RPhdf3U+7+iYMIG7wuv7qfd/QeF1/dT7v6GDCBu8Lr+6n3f0HhVf3U+7+gMIG7wqv7qfd/QeFV/dT7v6AwgbvCq/up939B4VX91Pu/oDCBu8Kr+6n3f0HhVf3U+7+gMIG7wqv7qfd/QeFV/dT7v6AwgbvCq/up939B4VX91Pu/oDCBu8Kr+6n3f0HhVf3U+7+gMIG7wqv7qfd/QeFV/dT7v6AwgbvCq/up939E+E1/dT7v6AwAb/Ca/up939B4TX91Pu/oDABv8Jr+6n3f0HhNf3U+7+gMAG/wmv7qfd/QeE1/dT7v6AwAb/Ca/up939B4TX99Pu/oDABv8Jr++n3f0HhNf30+7+gMAG/wnUe+n3f0HhNf30+7+gMAG/wAJr++n3f0HhNf30+7+gMAG/wAJr+6n3f0HhNf3U+7+gMAG7wqv7qfd/QeFV/dT7v6AwgbvC6/up939EeGV/dT7v6AxAbPDa3up939B4dW90O7+gMYGvw+t7od2Q9BVX7od2BlA0/g1fdDuR+HU6x7gbUi6RCRdI0oSLJEpF0gISJSLJFkiiqRKiXSJsBTEmxfEnEClgsMxDEClgxGYhiAvEnEviGIFMQxGYhiAvEMRmIYkC8QxGYk2AViGIywWCqYhiMxDEgXiGIywWCF4hiMxJxAViGI3ELAKxJxGWDEBdgxGYk4hSsQxG4hiQKxJxGYhiELxDEZiTYKViGI3EMSBWIYjcQxAViGI3EMQFYkYjcQxAViRiNxDEBLiVcR7iVxCEuJRxHuJVoBDRRoe4lJRAS0LaHtC2gBIukQkXSNiUi6QJF0igSJSJSLJFEJFkiUi6iVVEicRqgXUBgQoE4M0KBbAYM27ZO7ZpwJUBgy4Bgat2Tu0MVl3bDds1bsN2MGXBhgzVuw3YwZcAwZpwDdomDNgw3bNW7sGHwBlxDE07sMCYM2IYmndoN38DEZsQxNO7fQN2+hBmxDEdUlTpLzzjH+WZJ7S0sG0pOf8IgdiGAultLSzXGeD/wByNcHTqLyTjK/RgJwBQNO7J3YVlwZOBp3ZO7QMZd2w3bNW7RO7IYy7sndmrdoN38EMZd2G7Ne7+AwCsm7DdmvAN38AZN2G7Ne7DdkRj3ZDpmzdlXTQGTBlcTW6ZV0+qAy4lWjQ6ZRxsEIaKuI5oq0Ahoo0OaKNBCJIW0PkhUkURFF4orEZFHQWSLpEJF0VUpF4oiKGxiVBGAyMC8YjYx6msUqMPgYoDEki6h1LikqBZQHqKJSGBKh8EqA/H0BQGBOBOA7AFEYE4Bu/gfgGIxSMPgMDRjcjAmBGAbsfgGPwMCMPgMDRgRiTDGfEMDRgGFiBGAYfBowIxARuzm67XqmnToq75OXQ07Y1f4unxi/PNcLehwtNpdRq6mFOMpSfE59dYf8ABMoylLKbbfyRulwf/bO2v6Z1eF3Vhl0uWj/TOs4XqQ7nH9z+2/8A1df04DpfBChOHmhJp/DPTR/puag3OpFNGKjsWtqKklDyxX7mJ3Evi6/orQ7WlBxhqnlDll6o7lNRqQUoNSi+TRx6mwtRFySjlb1XqV2Tq3o625qN7pv19Gb56lS83n67mBKgPUU+K5EqBsIw+AwNCgTgQxnwJw+B+BKiDGfd/AYfBpwDEi4zYBgacfgMSGM6pk7s0YfBGBDGZw+Crh8GrG5DgExkdMo4GxwKSjwBjE6YuUPg2ygKlGxEYZwFNGycBE4lRnaKNDWikkVCJIXJD5CpIqKRGRRSIyJ0Fki6RVDIriaVeETRCNkUhE0QRqCYRGqARQ2KNyKhRsXUbkqIxIuKoodS6iWSJSArYLDMScRgXiGI3ELDFKxJxG4hYgViGI2wYkCsQxG2DEYFYInEZiGPwTAvEMUMxDEgViTiNxM20Zuhs+vUTs4wdiDiOHiO1puXGlTdv7JnpaFKnTglTgor4R57YyjHTx9zd5PqejpNOCPneXq2vV4JJF7k5Mgq+Zzehdu6syElbgrAkBERax5T+oNEqGqVSCtCpxX8nrDjf1LjLRRX7lI347nTn5ZvI2FW/I0CUn5qbxZ0sbHB/pRuVTUx+Fw7npMT3T48c+E4k4jcSbBSsfgnEZiTiRSsScRmJOKIFYhiMxJxMhWAYjVELAJxKuI/EholRncSriaHEXKJEZ5REyjc1tCZxAxziZ5xNtSPAz1I8CssUlZi5IfUQmRqMkyFyQ6QqQCojIlIjEdkXiOpoVEfBGop8FwHwQqmjRFXNxV4IbFFYjYo1IqYoukCRZI0BIvYlIskXFVsTYtYMQIsRYuo2LWJgXiTiXAiqYhYvYLEFLE2LWCxBWwWLWJsBSwWL2CxBSxk2tRdbZmohHm4XN1iHFSTT5PgSweY2JFypRl8ncr67TaOC3tRJ9FxZxHCeloShQlbzNXXS5yq0nLJq87c3LifPvH6rpPJ+Jkele39JJ2i2zVT19KpByi7pK54mk3JJunFX6HrNh6dfg5SXGROuJzHXx+Trq4pqdvUqE7Rg59StH+ptPUmoypTivVnG2hSen1MqaXrzM0pVaeoVNKM7/Bqcc2OfXl7le5oailqaedGakv8o4v9TXWnhL0ysznaGTjXU6TcJx5rqdrbNPfbGqTa81lI5/n89Rv9/visP9I0/PqaluHBX7nprHF/pOEVs6bv5pVHfsju2PZPjhz8UsFhiQWCqWJsWSJxIqlibF7BiQUsFi+IWMilgxL2sBAuxVxHWKtEQllZIbKJRoIRKIqSNMhMlcIyzXMzVEbKiM1RFSsVRGeSNdRcTLIrJUhUh0hUioVEYhcRiOyGRNFPmZ48zTTRqK0U0aIIRT5miC4nSKbFDYoXFDoo1FSkXiiIjEaEpEpAiyRVCRKRKJCosBJNgK2JJAggLE2JMitgsWsFkBUmxNibEFbBYtYCCLfAvUQc9PUiubi7DgsSxXn9maZanQqEuDxsKnslR8tnFr1XqbtNLcV6ySss3w/ub1NSjdcj5lt56r0ccTqPM0dkVKtfhkqa5ykei0tONGhu1yiRVq4JL1b4IvFq3PmZ66vTrx45ywbQ2XHWf6kXjUX+THpdlwztVlOEl6ejOzOrGko5uybtcbF+ondkw68ct1lhs2hGpC1K1l+pCtutU9nVY9VY6WZzdsUfyNO6d7XaMy7ZrHXOc3E6eE1PR0qNqdK2UlH1sjrpGHZGldDSwU3eSXY6Fj2+Kf4uXf0WIsWA2wrYmxJJBFgsTYmxFVsFi1gM0RYixYCChDLNEEMUaFyQ5i5EZJkhUkPaEzXMIzzRmqmqojPVXArNYqqMs+ZsqrgZKnM0yTIVIbIVIIVEZEXEZE7IbDmaaZmhzNVPmbitNMfART5GiB0inQQ2IqA2JpTEWRWJdGhZFkQiyKoJBEhQAEgQSQSZqAlEAiKsQACiQIJIAAAipAAIONUk47TrU2uEvNcfvVTgxW04ulrYVVynG390IqVU1e/A+d5ec7rv4+sheo31WlU1EZY48InHWprqzlUd/g319eqkXSTUYLm2Y4x00+O/tbnePMczJ7Z7v6vp0dnx1Ou09SlXq5Q9OHE3bNrVEnp6/wCuDtfqczTa2lRu9PUvZcYy4M209VDVSjUjZSsY6ldOOv8At10Zas3LU04+mQyFZOC4ldBHeaqc3ZqKsv5M8cfrrGu+sjoU44xSLogk98n5mPPbt0AABEgAEEklUSZVIABFFiCQMipDLFWBBSXIuUkZZpUhM0OkKnyDLPNGeryNMzNULGax1eRkqczZV5GSoaZJkKkNkKkVCYjYiojYnVDYI1UjNTNNI3FaaZogZ6RogdIp0eYyIuPMZE2pqLIquRdGlWSLEIkKlEgAQEkEgBBJBFSAAREpgQWCoJIAyJAAIAkgHwJVY9q0d9opNO0qfnX9jzb1N6bT5mjbW25VdatDppWp8pyX7vg5i4OzPJ5bLTT5bMpzpRlUqyWXqnyLw2dRSSWsgrdQp06langnw6FY7IrzlwS7nLf/AK6TP60+nsXSzqpvUSm/9rsUnRjoaripOSXI2UNn1dJFyb4/Bj1cs/13uZl2tdSSfMaY61qm38HX2FNz0cpvnKbPMU06klH0PV7Iio6O3ydfFM6c7bW8CESj0USAAZAiSCQAlEAjNFgQARUkABlUEMlkMgqyki7KyIzSZC58hrFTCM8+RmqGmfqZqgZrJV5My1DVV5My1DTBEhUhshUioTEbEVEbE6odTNNIy0+ZrpnSK0UjREz0jRE3FOjzGxFQGRNxTVyLoWi6NC6LFYlitJAAQRIABBIEEhUEkAZRIJgCIqQABRICNRq6GmjerUUfj1OTqdvOzWnp2/3SOfXcn1NdxtJXbsjz+09p1KspU6Tcaa6c5GWOr1GplKVWtJpLkZ692nY4d+TfUS1zNNDPaTm/Tkba8GmxGmWGsX8o6Wqp3WSOFvtqT/Fio6qVJ8HxHratWMrp8DJOmmJdL5GQnVjqS2tUqRs3/HEyyqyrSu2Z400aaUfMrDJF23626Olxv1NWk11TS7VlBv8A0mknH/5I08cUrmGpUVTW1J+iajczzf8ALWuvXL2id0muTJOTsvXcFQqu6/bI6yZ6+bqSpJIJCglEAQSAAQSuRJVMsZUIkgCKCrJIZKijKyLSKMiKMTMcxMghEzPVNEzNU9QzWWr6mSqa6nIyVSsESFSGyFSKhMRkRcRkTqh1PmaqZkhzNdM6RWimaIGekaIG4p8RkRUeY2JtTIl0LiXRoMRYoiyKqwEEhQSQSECJIAgAJIIqSs5xhG85KK6tla1WFCk6k3ZI83rtZU1VV5PGMeUfQ5d9zlLcdnU7X09FWg95LpHkcnVbW1NdWT3cekTFLhG9hlDT1K7/ANOEmlxbSbsebry9VnaTKTd3JttlZehr1tClRmo0pTkrcXJevqZJc1Y5o10Xu6Cs+L9BFV3dhOq1a0lCLteVuCuc57WrtNqlD/JcNdCcHCqpG+M1Ugnc42n2hKteNaKT9GlY6dGDUc48YmOo6cX+FKlKzfARKB1XSzjx4iKmjf7USVu8sUaY+jTcXdmijo2n5hlRKmuBLVnKtbU7nTtr9T4IToaUdx5/1SlfmZ5Xq1kn2NemcE1HOOXrHJFkxjrrTZRx/SndM0aTatXSSdOadSHOzfFAqE6y8lOT/gyVtNqI1pN0Kij1xfE1OsZ9vT6TWUdXDKlLj6p80aTxtOpOjVU4SlFrodzQ7YhNRhqHafLL0O07361OnWAhO6unwJK0CSCQAlEAjIsQBFyKCGDIZkVZRlmUZEUfNi5l5Cpv0CFTZmqGibM1RlZrLUMlTmaqnIy1OYYJkKkMkLkVCUMiKiNidUNhzNUOdzJE1QZuK00zRAz0x8TpFPiNixMRkTanIuhSLxNBiZZC0WQUwCqZYqpJKgBYCLgRFiAObtjWbmjuoPzzXZGeupzNowbW1v5FXCL/ANOHL5Zz1PzSTfB8iJPh/YpB+dfweDq/q7WF01Zne2HXjPSVKEoryyvaK4yv1OVoNDV1knZONK/Gb5HejU0uz6apQ4L1fq/5OPk6/iO3i4tu34XqdmR1PHLdv0VrnK12y3paaqb1Ts+KtY7ctdp6dKVWVWOK+TFOjPatDfZOmv8Ay4v1/k589Wffjt34+c9T28zqNLGs1Ju9vQTLR+RpWOjOLjJxkrNcGhbSPTrxY50aDhO7R2dk6mlaVGtNJvjHIyOHDiZdRSeSmlwX6l8Cz9TF5v5uvVwoqL+B6pRkjj7F1rhKOl1F0pf8uTPQKKXoeXvebj3eOzubGR0km3Yx66VDTU8q0lFtcF6s36ytDSUZVJWv+1dWeQ1lSVau5zd5S5m/Hz+vbHl7nHqDU7QjF/8AD0vPb9UuP+A/p/cvaVWvrae+lGGSUldZX6FIUL8WdHZNOFDU5qDnUf6UduvUebj31HptHWrVFeVDdwfXg3/Y1Sk0r24GOT1lOKqyjGol+qEXxRqoV6deiqlN3TPHZ/L3xn1mghqoOUVGNW3CXX+TgzjKEnCosZR4NM9Pxg72vHp0M+v0i1VO9O28XJ9fg6cd56rl5fHvuOTpNfqNJwUsoL9suR3dFtGlq/KvLUtdxZ5maWbU/La6t8lYVHSaqRbi8uFmemdWPNLY9mBz9nbRWoSp1WlU/wCp0DpLrpLqQIIuBYgi4XJRJVsGVZkQykmWbFyZBWQqYxsTJkQqoZqj4D6jM1RlYrPVMs+ZoqczNPmVkqQqQyQuQQmIyIqIyJ1Q2Jph6GWJpg+BuK1QY+BmpsfFnSK0RY2LERY1M2pyZZMWmXTKGJlkxaZZMovclMpcm4Uy4XKXC4F7k3KXC4EzmqcJTk7KKuzyuq1D1FadV+r4fwdnbWow0qpp2dR/4PPX426nl83W3Gequ2ny4Kx0dkbOVear11/pw5J/u/8AozbOoxrajKp/y4LiurNWr2jUpSkqa8voePu35HTx8z/atO09orTpU6dkrckcDVa+pWdru3RDI0a+um5Pixj2fKjHKcbLqTmTlvq9dM2jpOrUiqrtDnY7L18tM4wptS9FE413fym3RUlGe+rvlyHU36c9WfEaqM415uorSk7sRfhw5mjV197NX5R5OxnfH4Nz449fQvW6vcj/AKk8vUrJ+qVzTJUqUpVZVMm3H1PU7M1S1WjjOUllHhPj/k87FxqS8tk+TVys606EalGlUsqn6jHfP6mOvj8n4utW1tctTXeDvShwj8/Jz408nlIilTeeTbZpivixqTJkc+ur1dqsYWskr/Jpo6qOld0vN1EpJHS2Y9JKWFalCc37lcz1fTfjm9OrodbDU0k01l0LSoOFR1tOkpP9cPSXz/Jlr7MVG9fQeSXN0lyf8dDRs/WR1NLimprg01Zpnnv9x7Jb8rTGTau1b4JT6AznT17jqNykrp8TMmt7J9Z9sRp7+EopZv8AV8nNq2goxfPnY27TlGVW7OdU/Wk+h6ePjweT/amRqOju3Tk1K97nq9HqFqNNCp6tcV0Z47L/AFUr8lY7Ww9TjVnRb/VxX8nXmpzfbvXIK5BkbdFrhcpkDZETchsi5DZFDYuTJkyjZERJiZMvJipMiFVGZqjHTZnqMM1nqPiZ5Dqj4CJFZpchUhkhcghMRkRURiOqGRH05eghDKb4morXTZogzLFj4y4HSDTF8RkWIixkZGpVPTLpiVIumaU1MtcUpFlIoYnYtcUmTcBiZNxeROQF7k3F3IyIrh7YrZ63H0grHOnFq6QzUyc6tST5uTCit7KKfNP/AAeHq7dY+1t06lRopS/VJXYqvXjbFxvYbq53grPkcypL0ucZ79u19enX2dtOlp4veRsuonV6ye0qyjTTwXJHJm3K0b2SN2ztfS0N24ZzfK/oLz/LU72fmte4WlgnNLL5F0F+RVcpO1OP+Qp6ha3UXqvg+ZfWThCG6pNJEW59Ir1FVrSlbhyQr0tbiRfgQn8nR57dqWQ7dCPgpNvkvXmVFZzW8iox81+ZfdVZyznwhFWj8hTpqLTtdv1HuSUeDsQUjFIvdJ2VmhUpWZLkrc+LAb6cEXjaMoyg7SXIRTl5XdJr5L3ptckv4YV39m656i6k7SXodBRg5ZYrLqeRzdN50pSUl1ZdbX1K4RnxOPXj9+nq483r27u0tRW02Lp2afU506ueq3jtfqZ46qrXx308rDnDpxE5w66/SurlvJRkZarWa9OBeut20r8zHWk74x5vgdeZ6ebu+0wleUpLhxsjZpqu4qRmucTFTSvw5JD1w9eJpl7CM1OKkuTVybmHZtZz0NPquBqyOjrplyMijkRkRV8iHIpkQ5EFmykmVcijkETJiZyJlITJkS1WbM9RjJyM85BkqoxMi83xFSDKkhci8hcghMRiFIvE6hsWMixSZdFg0wfAdCRkjKw+LTNxWqMhsZGWMxsZGtVoUhikZlIupGtGhSJUhCkWyLoepE5CMicgp2QZCVInIaHZCtTV3enqSvZqPAjIwbUq3hCmut2Y66yajm343t/Jq01NQi5v9y4fwZru3wPrVFuVCL4o8PTXE/lnrVXKbS5CKjjF3ZN2m2mjLUk3Jt8WItqXNv4KMi9y8KTlxtw+Ssr0JVHJYya+TYrt35vqxcI4qxZsJavzIcn0RW9lxYf5KiW+IJWVyv8AIOf/AEIJcnfgRf0SX9il7u9yMr8ih100rlOb9SqdgcunMCyvHhz+SynK3GP+SnNceZCdndAOUsk1a1xTgnLgS7Lj1KZ2kStQ6nNwkdDTVM7cePJnMTUmjZpLxqrozFdOTNpeRRn14HOybeT5yOptOm5aGTXo0zkU1L04/P0a4+MeSZTlPHyx4y9fgvBvrfqysKSS4/8A6XtZcVZG3N3tkT/4S3STNuZy9lz/AOGf/qNmZqOkvpocyMxOZDnxC6dmQ5CcyrmDTnMW5FHPiLcyIvKQqUiJSFSmRBOQiciZSFSYRVsXJlmyjZUUkLkXkxUmEJTGRYlMYmdA1MumKTLpgNTLxlYSmWTNDUpoZGZjUhkZ2NarWpl1J9TIqiLqoXRqU2WzMqqE7xDRqzJz+TLvETvC6rSp/JORmzuG8Jo05h+FHVeabcZW4WEQlnNROhSlZroeX/yPJf8AWO3j5l91x9doZaOcLzUoy/sznVXJu0TftPUOvqZXd4x4I58pWRym/wAs9ZL6Kl5Y2vxKJJlZviTCLZtk2MIrixiZSVopIi9uBYx19NuTlfgJytyZKl1CG34Blf6F5WfEG+gVdy/sVb6lXLiRYCX/AHAhsi4RZMlMWnYnLiAxyKyb5kZEN9Shil/2yjWTK3foWg/MiLDafA6mndopxjxOdTcVLzHSo6ylDy2Xwc+nbk2tNTpyptcJrocuPl8r9OB0d7k3a1jDqYqFd9JK44/o8nv2E4ri7t9CLNu8iqcfktf/AGnRxdDZ8rUWvk17w5ukqWUlyNG8NxWnMMzNvA3gVocyrmIc/krvANDmUcxDqFXUIhspi5TFuVyrYFmyjZDZVsIGxbZLZRsIiTFSZaTFyYCUy6YlMYmbDUy6YpMsmUNTLJikyyYDUybi0yci6G3JUhWRNwpqmy28E3C5dD94G8EXJyGofvAzEZEx80kiarqaGF7z6jNbV3FHnZvghMdVChBL16IyaqpOu85cF6LoeK711terZzzkZJvqZqjNEzNN8TpHClpNsfBJCY3uMV1Fvoiopvc61RekeBVyM9Jyp1ZxmrN8S8maYMjJ9i+Qi/qWUgHX9Cbr0E5F4sC+XQm91a4vNL0Ji7gMXFENAv4BtEFWRyJkUKL3XMq5FW+JVt+lgL5di0JeZCW30CMmBvmiIPjxZbmkLfBmHR0NNO8cWGtheMZpfpdhFCXmRsi23wXDqY+XXT7MZ6WnrT/5dKpL+IsetBq7X/Hml88D0ez6290kHazSsxlaaUWZvlv9NzwTN15Z0qlGXnjjf5DI17TfJ/Jz7nfi7Nefufm4bkGQq4XNsm5EZC8iMiBmRDkUuRcC9yrZVyIbAs5FWyrZVyCJbKSZDkUbAJMW2S2UbASmXTOf+VPpEn8up0ia1HRTLpnM/NqdI9ifzqvth2Y0dRMlM5X59X2w7MnxCt7YdmNV1kybnI8Rre2HZ/ZPiNb2w7P7Lo69ybnH8Sre2HZ/YeJVvbT7P7GjsZE3OP4nW9tPs/sPE63tp9n9jR2MgyOP4nX9tPs/sPE6/tp9n9jR2chlK9+HM4Xidf20+z+xkNs6iCsoUuz+zPV2ZF5uV6CEFF5S4v5K1pp8jhPbeqf7aXZ/ZR7X1D/bT7P7OP5rp+46tWRm5swS2lWlzjDs/sqtfVvfGHZmsYvUdOKHQXU5HiNb2w7P7J8Urpfpp9n9lk9mw/U8dc7crK4uc7vhyMr1NRylKyvJ3ZXfS6I1WWu/Qm5k/In0iT+TPpEg2XGQZg/Kn0iT+ZU6R7BHRfwy0X6HN/Nqe2HZk/n1fbDswOjOXDgyinZ2ZgetqP8AbDsR+ZU9sOwV0nJMo2v7GD8up0iR+XU6R7BG5yKN8eBj/Jn0iH5M+kQrXk+pMZceZj/Jn0iH5M+kQO7B3pxfwVkcqG060IqKjCy6p/ZL2nWf7afZ/ZMa/UdeizpUOMUzy0dp1o8o0+z+x0NuaqHKNLs/szea3O5Hstm6jdznTfrxRsq1bp3PCf8AiDVqakoUU18P7Gv+qNc+dOh/7X9nO+KuvPm5kx6LX+am/g5eRzZ/1Fq5xadOhx/2v7M3itf20+z+zrxLJlcfJ1OrsdvILnE8Vr+2n2f2Hitf20+z+zbnrt5BkcTxWv7afZ/YeK1/bT7P7BrtZBkcXxWv7afZ/YeKV/bT7P7A7LkVcjj+KV/bT7P7DxOt7afZ/YR1nIq5HK8Sre2HZ/ZHiNb2w7P7A6jZVs5viFX2w7Mj8+r7YdmB0GyjZhetqv0j2I/LqdIgZwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//Z",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"560\"\n",
       "            height=\"315\"\n",
       "            src=\"https://www.youtube.com/embed/__cOYPifDWk\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x75c691db97d0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('__cOYPifDWk', width=560, height=315)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Table of contents**    \n",
    "1. Import required packages and set up some stuff    \n",
    "2. Retrieve First-Level results    \n",
    "3. Displaying subject Effects-Of-Interest z-maps    \n",
    "4. Specify the second-level model    \n",
    "4.1. Design matrix for 3x3 repeated-measures ANOVA   \n",
    "4.2. Model specification and fit Contrasts  \n",
    "4.3. Contrasts   \n",
    "5. Computing contrasts and plotting result maps    \n",
    "5.1. Faces > Scrambled \n",
    "5.2. Famous > Unfamiliar  \n",
    "5.3. Stimulus x Immediate Repetition interaction  \n",
    "5.4. FWE correction using non-parametric permutation testing    \n",
    "6. Summary results  \n",
    "6.1. Using atlasreader package    \n",
    "6.2. Nilearn's report   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "## Import required packages and set up some stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The conda environment used for this tutorial is available here: https://github.com/MRC-CBU/COGNESTIC/blob/main/mri_environment.yml \n",
    "\n",
    "import os\n",
    "import glob # to search for files using regex\n",
    "\n",
    "import pandas as pd # for data manipulation\n",
    "import numpy as np # for numerical operations\n",
    "\n",
    "from bids.layout import BIDSLayout # to fetch data from BIDS-compliant datasets\n",
    "\n",
    "import matplotlib.pyplot as plt # for basic plotting\n",
    "\n",
    "import nibabel as nib # NiBabel, to read and write neuroimaging data, https://nipy.org/nibabel/\n",
    "\n",
    "# Nilearn modules, for the analysis of brain volumes, plotting, etc., https://nilearn.github.io/\n",
    "from nilearn.plotting import plot_glass_brain, plot_design_matrix, plot_contrast_matrix, plot_stat_map, view_img, view_img_on_surf\n",
    "from nilearn.glm.second_level import SecondLevelModel\n",
    "from nilearn.glm.thresholding import threshold_stats_img\n",
    "from nilearn.datasets import load_mni152_template\n",
    "from nilearn.glm.second_level import non_parametric_inference\n",
    "\n",
    "from atlasreader import create_output # For generating result tables https://github.com/miykael/atlasreader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MNI152 template will be used as a backgound for plotting results\n",
    "mni152_template = load_mni152_template()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve First-Level results\n",
    "\n",
    "For the group analysis, we will use the single-condition contrast estimate (beta estimate) maps for all nine conditions. Because we saved the results in BIDS format, we can us PyBIDS to retrieve the subject-level results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Set up the paths to the data and results folders\n",
    "fmri_data_dir = '../example_data/FaceRecognition/data' # data in BIDS format\n",
    "fmri_results_dir = '../example_data/FaceRecognition/results' # results in BIDS format\n",
    "\n",
    "fmri_group_dir = os.path.join(fmri_results_dir, 'group-level') # where group results will go\n",
    "os.makedirs(fmri_group_dir, exist_ok=True)\n",
    "\n",
    "# --- Set up the BIDS layout\n",
    "layout = BIDSLayout(fmri_data_dir, derivatives = True)\n",
    "\n",
    "# Attach the results folder to the layout. It must complay with BIDS standards. \n",
    "# And must include dataset_description.json file!\n",
    "layout.add_derivatives(os.path.join(fmri_results_dir, \"first-level\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Displaying subject Effects-Of-Interest z-maps\n",
    "\n",
    "To check how the first-level results look overall, it is helpful to display effects-of-interest for all subjects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eoi_maps = layout.get(desc='EffectsOfInterest', extension='.nii.gz', return_type='file')\n",
    "print(*eoi_maps, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects = layout.get_subjects()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=4, ncols=4, figsize=(14, 14))\n",
    "\n",
    "for i, stat_map in enumerate(eoi_maps):\n",
    "    plot_glass_brain(stat_map, \n",
    "                              title = 'sub-' + subjects[i],\n",
    "                              axes = axes[int(i / 4), int(i % 4)],\n",
    "                              plot_abs = False, \n",
    "                              display_mode='x')\n",
    "fig.suptitle('Effects of interest' + ', unthresholded z-maps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify the second-level model\n",
    "\n",
    "At the group-level analysis, we also use a GLM. The outcome variable is the beta/contrast estimate from each subject, and the predictor variables typically include group-level factors such as experimental conditions, subject-specific regressors (in repeated-measures designs), group-specific regressors (in between-subject designs), and subject-specific covariates (e.g., age, gender, or behavioural scores)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Design matrix for 3x3 repeated-measures ANOVA\n",
    "\n",
    "In this example, the design matrix we generate will incorporate both **within-subject (conditions)** and **between-subject (subjects)** factors. Each row in the design matrix corresponds to a specific observation, which in this case is a beta estimate from a given condition and subject, while each column represents a predictor variable.\n",
    "\n",
    "The number of rows in the design matrix must match the number of first-level result files that will be entered into the second-level model. The order of the rows in the design matrix must match the order of the provided files.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify which conditions to include in the analysis and retrieve their effect files from the first-level results.\n",
    "conditions = ['IniFF', 'ImmFF', 'DelFF', 'IniUF', 'ImmUF', 'DelUF', 'IniSF', 'ImmSF', 'DelSF']\n",
    "\n",
    "effect_files = layout.get(desc=conditions, suffix='effect', extension='.nii.gz', return_type='filename')\n",
    "\n",
    "# print to see if it found what we expexted\n",
    "print(f\"Found {len(effect_files)} effect files:\")\n",
    "print(*effect_files[0:10], sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note we have reordered conditions so that the \"stimulus\" factor (Famous Face, Unfamiliar Face, Scrambled Face) rotates slowest, and \"presentation\" factor (Initial, Immediate Repeat, Delayed Repeat) rotates fastest in this 3x3 design."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Create an empty DataFrame\n",
    "df = pd.DataFrame(columns=conditions + subjects)\n",
    "\n",
    "# Populate the DataFrame with 0s and 1s\n",
    "for i, condition in enumerate(conditions):\n",
    "    # Filter files based on condition\n",
    "    condition_files = [1 if condition in file else 0 for file in effect_files]\n",
    "    # Add a column for the condition\n",
    "    df[condition] = condition_files\n",
    "\n",
    "# Populate the DataFrame with 0s and 1s for subjects\n",
    "for i, subject in enumerate(subjects):\n",
    "    # Filter files based on subject\n",
    "    subject_files = [1 if f\"sub-{subject}\" in file else 0 for file in effect_files]\n",
    "    # Add a column for the subject\n",
    "    df[subject] = subject_files\n",
    "\n",
    "# Display the design matrix\n",
    "design_matrix = df\n",
    "ax = plot_design_matrix(design_matrix)\n",
    "ax.set_title(\"Second level design matrix\", fontsize=12)\n",
    "ax.set_ylabel(\"stat maps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first 9 columns capture the mean of each condition, while the remaining 16 columns capture the mean of each subject (which we do not care about, and removing this variance can improve statistics for contrasts over the first 9 columns by reducing the residual error). \n",
    "\n",
    "### Model specification and fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify\n",
    "second_level_model = SecondLevelModel() \n",
    "# fit\n",
    "second_level_model = second_level_model.fit(\n",
    "  effect_files, \n",
    "  design_matrix = design_matrix\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contrasts\n",
    "\n",
    "We can specify some contrasts that might be of interest. The first two are T-contrasts (one row), the last is an F-contrast (more than one row):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_columns = design_matrix.shape[1]\n",
    "\n",
    "contrasts = {\n",
    "  'FacesScrambled':   np.pad([1/6,  1/6,  1/6,  1/6,  1/6,  1/6, -1/3, -1/3, -1/3], (0, n_columns - 9), 'constant'),\n",
    "  'FamousUnfamiliar': np.pad([1/3,  1/3,  1/3, -1/3, -1/3, -1/3,    0,    0,    0], (0, n_columns - 9), 'constant'),\n",
    "  'Stimulus X Immediate Repetition': np.array([np.pad(np.kron([1/2, 1/2, -1], [1, -1, 0]), (0, n_columns - 9), 'constant'),\n",
    "                                     np.pad(np.kron([             1,   -1,   0], [1, -1, 0]), (0, n_columns - 9), 'constant')])\n",
    "}\n",
    "\n",
    "for contrast_id, contrast_val in contrasts.items():\n",
    "    plot_contrast_matrix(contrast_val, design_matrix=design_matrix)\n",
    "    plt.suptitle(contrast_id)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there is no particular scientific reason for the specific F-contrast interaction above, other than to demonstrate how to construct such contrasts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing contrasts and plotting result maps\n",
    "\n",
    "We could then loop through every contrast like below, using the same threshold:\n",
    "```python\n",
    "for contrast_id, contrast_val in contrasts.items():\n",
    "    print(f\"\\tcontrast id: {contrast_id}\")\n",
    "    z_map = second_level_model.compute_contrast(contrast_val, output_type=\"z_score\")\n",
    "\n",
    "    thresholded_map, threshold = threshold_stats_img(\n",
    "        z_map, \n",
    "        alpha = 0.001, \n",
    "        height_control = None, \n",
    "        two_sided = True)\n",
    "    \n",
    "    plot_stat_map(\n",
    "        thresholded_map, \n",
    "        threshold = threshold, \n",
    "        bg_img = mni152_template,\n",
    "        display_mode=\"ortho\",\n",
    "        black_bg = True,\n",
    "        cmap = 'hot', \n",
    "        title = contrast_id)  \n",
    "    \n",
    "    fig = plt.gcf()\n",
    "    fig.set_size_inches(10,3)\n",
    "    plt.show()\n",
    "```\n",
    "but let's look at each contrast separately. \n",
    "\n",
    "### Faces > Scrambled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the z-maps for the contrast\n",
    "z_map = second_level_model.compute_contrast(\n",
    "  contrasts['FacesScrambled'], \n",
    "  output_type='z_score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholded_map_bonf, threshold_bonf = threshold_stats_img(\n",
    "  z_map, \n",
    "  alpha=0.05, \n",
    "  height_control='bonferroni', \n",
    "  two_sided=False)\n",
    "\n",
    "print('Bonferroni p<.05 threshold: %.3f' % threshold_bonf)\n",
    "\n",
    "plot_stat_map(\n",
    "    thresholded_map_bonf, \n",
    "    bg_img = mni152_template,\n",
    "    threshold=threshold_bonf,   \n",
    "    display_mode = 'ortho',\n",
    "    black_bg = True,\n",
    "    cmap = 'hot',    \n",
    "    title = 'Faces > Scrambled  (Bonf. p<.05)'\n",
    "    )\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(10,3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to one-sample t-test, we have many more degrees of freedom in this 3x3 ANOVA model, so more power. However, this power comes at the price of assuming that the error is spherical, i.e., that there is equal variance across the 9 conditions, and equal covariance. If this is not true, the statistics are biased, and hence the one-sample t-test approach is safer; see [Rik's Stats](https://github.com/RikHenson/CognesticStats) for more details. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interactive plotting\n",
    "plot = view_img(\n",
    "  thresholded_map_bonf, \n",
    "  bg_img=mni152_template, \n",
    "  threshold=threshold_bonf, \n",
    "  colorbar=True, \n",
    "  title='Faces > Scrambled  (Bonf. p<.05)'\n",
    "  )\n",
    "\n",
    "plot.open_in_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Famous > Unfamiliar\n",
    "\n",
    "Next, what about the famous vs unfamiliar contrast. This time we will use a one-tailed FDR correction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_map = second_level_model.compute_contrast(contrasts['FamousUnfamiliar'], output_type=\"z_score\")\n",
    "thresholded_map, threshold = threshold_stats_img(\n",
    "        z_map, \n",
    "        alpha = 0.05, \n",
    "        height_control = 'fdr', \n",
    "        two_sided = False)\n",
    "plot_stat_map(\n",
    "        thresholded_map, \n",
    "        threshold = threshold, \n",
    "        bg_img = mni152_template,\n",
    "        display_mode=\"ortho\",\n",
    "        black_bg = True,\n",
    "        cmap = 'hot', \n",
    "        title = 'Famous > Unfamiliar',\n",
    "    )       \n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(10,3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see some medial prefrontal regions that are more active for famous faces (that participants were likely to recognise, and may trigger semantic/emotional responses).\n",
    "\n",
    "### Stimulus x Immediate Repetition interaction\n",
    "\n",
    "Finally, what about the interaction F-contrast we created, which asks whether the effect of repetition depends on whether the repeated stimulus is a face or scrambled face:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_map = second_level_model.compute_contrast(contrasts['Stimulus X Immediate Repetition'], output_type=\"z_score\") # Don't need stat_type=\"F\" option?\n",
    "thresholded_map, threshold = threshold_stats_img(\n",
    "        z_map, \n",
    "        alpha = 0.001, \n",
    "        height_control = None,\n",
    "        two_sided = True)\n",
    "plot_stat_map(\n",
    "        thresholded_map, \n",
    "        threshold = threshold, \n",
    "        bg_img = mni152_template,\n",
    "        display_mode=\"ortho\",\n",
    "        black_bg = True,\n",
    "        cmap = 'hot', \n",
    "        title = 'Stimulus X Immediate Repetition',\n",
    "    )       \n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(10,3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is nothing particularly convincing for this last contrast: unlikely that anything would survive correction for multiple comparisons, but some might be reliable when using more focused Region-of-Interest (ROI) analyses as in the [next Notebook](../notebooks/nb06_ROI_analysis.ipynb) (where we will also be able to plot the 9 conditions in order to understand interactions like this). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FWE correction using non-parametric permutation testing\n",
    "\n",
    "Nilearn's FWE correction using the Bonferroni approach (`height_control='bonferroni'`) is applied to the number of voxels. However, this method is not well-suited for fMRI data because neuroimaging data typically exhibit spatially correlated data points, which violate the Bonferroni assumption of independent tests.\n",
    "\n",
    "As an alternative, neuroscientists have developed **Random Field Theory** (RFT), which accounts for spatial correlations by applying multiple comparison corrections that consider the smoothness of the data (it can be thought of as a Bonferroni correction based on the number of **'resels'** - RESolution ELements - rather than the raw number of voxels). However, it makes various parametric assumptions about the data, including a minimal smoothness in terms of voxels, that can be violated. More practically, it is not implemented in Nilearn (to our knowledge).\n",
    "\n",
    "Another alternative, which is implemented in Nilearn, is to use non-parametric inference, specifically permutation testing. This makes fewer assumptions when correcting for multiple comparisons (though note that is assumes exchangeability, which can also be violated if the error is nonspherical under the null; see [Rik's Stats](https://github.com/RikHenson/CognesticStats)). The only downside is that permutation testing can take a long time to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_dict = non_parametric_inference(\n",
    "    effect_files,\n",
    "    design_matrix = design_matrix,\n",
    "    second_level_contrast = contrasts['FacesScrambled'],\n",
    "    n_perm = 100, # ideally at least 10000\n",
    "    two_sided_test = False,\n",
    "    n_jobs = -1, # Use all available cores\n",
    "    threshold = 0.001 # cluster level threshold; enables cluster-level inference\n",
    ")\n",
    "\n",
    "# Print the keys of the output dictionary\n",
    "print(out_dict.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is multiple images (maps), organised in a dictionary. \n",
    "* Voxel-level inference\n",
    "  * **t**: t-statistics\n",
    "  * **logp_max_t**: Negative log10 family-wise error rate-corrected p-values corrected based on the distribution of maximum t-statistics from permutations.\n",
    "* Cluster-level inference\n",
    "  * **size**: Cluster size values associated with the significance test \n",
    "  * **logp_max_size**: Negative log10 family-wise error rate-corrected p-values corrected based on the distribution of maximum cluster sizes from permutations.\n",
    "  * **mass**: Cluster mass values associated with the significance test \n",
    "  * **logp_max_mass**: Negative log10 family-wise error rate-corrected p-values corrected based on the distribution of maximum cluster masses from permutations. \n",
    "\n",
    "**We will focus only on the voxel-level inference.**\n",
    "\n",
    "To report the FWE-corrected maps, we could display the *logp_max_t*; however, these values can be difficult to interpret if you're not familiar with them. It might be better to plot and report a t-map, masked to exclude the voxels that did not survive the FWE correction.\n",
    "\n",
    "Let's create a new image displaying t-values for the voxels with a p-value < 0.05."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "masked = out_dict['logp_max_t'].get_fdata() > -np.log10(alpha)\n",
    "masked_t_map = out_dict['t'].get_fdata() * masked\n",
    "\n",
    "# save the masked t-map as a nifti image\n",
    "masked_t_map_img = nib.Nifti1Image(masked_t_map, out_dict['t'].affine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the smallest t-value that is above the threshold (for the colorbar; the maps themselves are thresholded already)\n",
    "threshold_fwe = masked_t_map[masked_t_map > 0].min()\n",
    "print('FWE (perm.) p<.05 threshold: %.3f' % threshold_fwe)\n",
    "\n",
    "plot_stat_map(\n",
    "    masked_t_map_img, \n",
    "    threshold = threshold_fwe,       \n",
    "    display_mode = 'ortho',\n",
    "    black_bg = True,\n",
    "    bg_img = mni152_template,\n",
    "    cmap = 'hot',\n",
    "    title = f\"Faces > Scrambled (FWE p <.{alpha})\")\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(10,3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we observe that the non-parametric FWE correction is slightly less conservative than the Bonferroni correction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary results\n",
    "\n",
    "### Using atlasreader package\n",
    "\n",
    "We can use ['atlasreader'](https://github.com/miykael/atlasreader) package to get summary results (peak table, cluster table, .png images of each cluster). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate and save atlasreader output\n",
    "outdir = os.path.join(fmri_results_dir, 'group-level', 'permutation', 'FacesScrambled')\n",
    "\n",
    "create_output(\n",
    "    masked_t_map_img, \n",
    "    cluster_extent = 20, \n",
    "    voxel_thresh = threshold_fwe,\n",
    "    direction = 'pos',\n",
    "    outdir = outdir\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display the peak table\n",
    "peaks = glob.glob(os.path.join(outdir, '*_peaks.csv'))\n",
    "display(pd.read_csv(peaks[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display the cluster table\n",
    "clusters = glob.glob(os.path.join(outdir, '*_clusters.csv'))\n",
    "display(pd.read_csv(clusters[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Some more plotting options**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the top 5 peaks' x values\n",
    "x = pd.read_csv(peaks[0])['peak_x'][:5]\n",
    "# sort the x values\n",
    "x = x.sort_values()\n",
    "\n",
    "# plot these peaks\n",
    "plot_stat_map(\n",
    "    masked_t_map_img, \n",
    "    threshold = threshold_fwe,       \n",
    "    display_mode = 'x',\n",
    "    cut_coords = x,\n",
    "    black_bg = True,\n",
    "    bg_img = mni152_template,\n",
    "    cmap = 'hot',\n",
    "    title = f\"Faces > Scrambled (FWE p <.{alpha})\")\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(15,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at a 3D brain using `plotly`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = view_img_on_surf(masked_t_map_img, \n",
    "    threshold = threshold_fwe)\n",
    "#view.open_in_browser()\n",
    "view.resize(1600, 800)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or use, for example, FSLeyes to plot and explore the result maps. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nilearn's report\n",
    "\n",
    "Nilearn has a built-in report generator that can create reports for all defined contrasts. However, a limitation is that it cannot generate reports for results obtained using non-parametric inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_level_report = second_level_model.generate_report(\n",
    "  contrasts, \n",
    "  title = \"Results of the second-level analysis\", \n",
    "  bg_img = mni152_template, \n",
    "  alpha = 0.001, \n",
    "  cluster_threshold = 20, \n",
    "  height_control = 'fpr', \n",
    "  min_distance = 8.0, \n",
    "  plot_type = 'slice', \n",
    "  display_mode = 'x', \n",
    "  report_dims = (1600, 800))\n",
    "\n",
    "second_level_report.open_in_browser()\n",
    "\n",
    "#second_level_report.save_as_html(file_name)"
   ]
  }
 ],
 "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.11.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "337px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "874.85px",
    "left": "2183px",
    "right": "20px",
    "top": "116px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}