{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Tutorial 01: Linear Models\nIn this tutorial you will set up your first linear model with BrainStat. \nFirst lets load some example data to play around with. We'll load age, IQ, and left\nhemispheric cortical thickness for a few subjects.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport nibabel as nib\nimport os\nimport brainstat\nfrom brainstat.tutorial.utils import fetch_tutorial_data\nfrom brainstat.context.utils import read_surface_gz\nfrom nilearn.datasets import fetch_surf_fsaverage\n\nbrainstat_dir = os.path.dirname(brainstat.__file__)\ndata_dir = os.path.join(brainstat_dir, \"tutorial\")\n\nn = 20\ntutorial_data = fetch_tutorial_data(n_subjects=n, data_dir=data_dir)\nage = tutorial_data[\"demographics\"][\"AGE\"].to_numpy()\niq = tutorial_data[\"demographics\"][\"IQ\"].to_numpy()\n\n# Reshape the thickness files such that left and right hemispheres are in the same row.\nfiles = np.reshape(np.array(tutorial_data[\"image_files\"]), (-1, 2))\n\n# We'll use only the left hemisphere in this tutorial.\nthickness = np.zeros((n, 10242))\nfor i in range(n):\n    thickness[i, :] = np.squeeze(nib.load(files[i, 0]).get_fdata())\nmask = np.all(thickness != 0, axis=0)\n\npial_left = read_surface_gz(fetch_surf_fsaverage()[\"pial_left\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we can create a BrainStat linear model by declaring these variables as\nterms. The term class requires two things: 1) an array or scalar, and 2) a\nvariable name for each column. Lastly, we can create the model by simply\nadding the terms together.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from brainstat.stats.terms import Term\n\nterm_intercept = Term(1, names=\"intercept\")\nterm_age = Term(age, \"age\")\nterm_iq = Term(iq, \"iq\")\nmodel = term_intercept + term_age + term_iq"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also add interaction effects to the model by multiplying terms.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model_interaction = term_intercept + term_age + term_iq + term_age * term_iq"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, lets imagine we have some cortical marker (e.g. cortical thickness) for\neach subject and we want to evaluate whether this marker changes with age\nwhilst correcting for effects of sex and age-sex interactions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from brainstat.stats.SLM import SLM\n\nslm = SLM(model_interaction, -age, surf=pial_left, correction=\"rft\", mask=mask)\nslm.fit(thickness)\nprint(slm.t.shape)  # These are the t-values of the model.\nprint(slm.P[\"pval\"][\"P\"])  # These are the random field theory derived p-values."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By default BrainStat uses a two-tailed test. If you want to get a one-tailed\ntest, simply specify it in the SLM model as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slm_two_tails = SLM(\n    model_interaction, -age, surf=pial_left, correction=\"rft\", two_tailed=False\n)\nslm_two_tails.fit(thickness)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, imagine that instead of using a fixed effects model, you would prefer a\nmixed effects model wherein handedness is a random variable. This is simple to\nset up. All you need to do is initialize the handedness term with the Random\nclass instead, all other code remains identical.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from brainstat.stats.terms import Random\n\nrandom_handedness = Random(tutorial_data[\"demographics\"][\"HAND\"], name_ran=\"Handedness\")\nrandom_identity = Random(1, name_ran=\"identity\")\nmodel_random = (\n    term_intercept\n    + term_age\n    + term_iq\n    + term_age * term_iq\n    + random_handedness\n    + random_identity\n)\nslm_random = SLM(model_random, -age, surf=pial_left, correction=\"fdr\", mask=mask)\nslm_random.fit(thickness)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "That concludes the basic usage of the BrainStat for statistical models.\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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}