{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Beta-Series Modeling for Task-Based Functional Connectivity and Decoding\n\nThis example shows how to run beta series :term:`GLM` models, which are a\ncommon modeling approach for a variety of analyses of task-based :term:`fMRI`\ndata with an event-related task design, including\n:term:`functional connectivity`, :term:`decoding `, and\nrepresentational similarity analysis.\n\nBeta series models fit trial-wise conditions, which allow users to create\n\"time series\" of these trial-wise maps, which can be substituted for the\ntypical time series used in resting-state functional connectivity analyses.\nGenerally, these models are most useful for event-related task designs,\nwhile other modeling approaches, such as psychophysiological interactions\n(PPIs), tend to perform better in block designs, depending on the type of\nanalysis.\nSee :footcite:t:`Cisler2014` for more information about this,\nin the context of functional connectivity analyses.\n\nTwo of the most well-known beta series modeling methods are\nLeast Squares- All (LSA) :footcite:p:`Rissman2004` and\nLeast Squares- Separate (LSS)\n:footcite:p:`Mumford2012,Turner2012`.\nIn LSA, a single :term:`GLM` is run, in which each trial of each condition of\ninterest is separated out into its own condition within the design matrix.\nIn LSS, each trial of each condition of interest has its own :term:`GLM`,\nin which the targeted trial receives its own column within the design matrix,\nbut everything else remains the same as the standard model.\nTrials are then looped across, and many GLMs are fitted,\nwith the :term:`Parameter Estimate` map extracted from each GLM to build the\nLSS beta series.\n\n.. topic:: Choosing the right model for your analysis\n\n We have chosen not to reproduce analyses systematically comparing beta\n series modeling approaches in nilearn's documentation;\n however, we do incorporate recommendations from the literature.\n Rather than taking these recommendations at face value, please refer back\n to the original publications and any potential updates to the literature,\n when possible.\n\n First, as mentioned above, according to :footcite:t:`Cisler2014`,\n beta series models are most appropriate for event-related task designs.\n For block designs, a PPI model is better suited- at least for\n functional connectivity analyses.\n\n According to :footcite:t:`Abdulrahman2016`,\n the decision between LSA and LSS should be based on three factors:\n inter-trial variability, scan noise, and stimulus onset timing.\n While :footcite:t:`Mumford2012` proposes LSS as a tool\n primarily for fast event-related designs (i.e., ones with short inter-trial\n intervals), :footcite:t:`Abdulrahman2016` finds, in simulations,\n that LSA performs better than LSS when trial variability is greater\n than scan noise, even in fast designs.\n\n.. include:: ../../../examples/masker_note.rst\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n\nfrom nilearn import image, plotting"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Prepare data and analysis parameters\nDownload data in :term:`BIDS` format and event information for one subject,\nand create a standard :class:`~nilearn.glm.first_level.FirstLevelModel`.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from nilearn.datasets import fetch_language_localizer_demo_dataset\nfrom nilearn.glm.first_level import FirstLevelModel, first_level_from_bids\n\ndata_dir, _ = fetch_language_localizer_demo_dataset()\n\nmodels, models_run_imgs, events_dfs, models_confounds = first_level_from_bids(\n data_dir,\n \"languagelocalizer\",\n img_filters=[(\"desc\", \"preproc\")],\n)\n\n# Grab the first subject's model, functional file, and events DataFrame\nstandard_glm = models[0]\nfmri_file = models_run_imgs[0][0]\nevents_df = events_dfs[0][0]\n\n# We will use first_level_from_bids's parameters for the other models\nglm_parameters = standard_glm.get_params()\n\n# We need to override one parameter (signal_scaling)\n# with the value of scaling_axis\nglm_parameters[\"signal_scaling\"] = standard_glm.scaling_axis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define the standard model\nHere, we create a basic :term:`GLM` for this run, which we can use to\nhighlight differences between the standard modeling approach and beta series\nmodels.\nWe will just use the one created by\n:func:`~nilearn.glm.first_level.first_level_from_bids`.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"standard_glm.fit(fmri_file, events_df)\n\n# The standard design matrix has one column for each condition, along with\n# columns for the confound regressors and drifts\nfig, ax = plt.subplots(figsize=(5, 10))\nplotting.plot_design_matrix(standard_glm.design_matrices_[0], ax=ax)\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define the LSA model\nWe will now create a least squares- all (LSA) model.\nThis involves a simple transformation, where each trial of interest receives\nits own unique trial type.\nIt's important to ensure that the original trial types can be inferred from\nthe updated trial-wise trial types, in order to collect the resulting\nbeta maps into condition-wise beta series.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Transform the DataFrame for LSA\nlsa_events_df = events_df.copy()\nconditions = lsa_events_df[\"trial_type\"].unique()\ncondition_counter = {c: 0 for c in conditions}\nfor i_trial, trial in lsa_events_df.iterrows():\n trial_condition = trial[\"trial_type\"]\n condition_counter[trial_condition] += 1\n # We use a unique delimiter here (``__``) that shouldn't be in the\n # original condition names\n trial_name = f\"{trial_condition}__{condition_counter[trial_condition]:03d}\"\n lsa_events_df.loc[i_trial, \"trial_type\"] = trial_name\n\nlsa_glm = FirstLevelModel(**glm_parameters)\nlsa_glm.fit(fmri_file, lsa_events_df)\n\nfig, ax = plt.subplots(figsize=(10, 10))\nplotting.plot_design_matrix(lsa_glm.design_matrices_[0], ax=ax)\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Aggregate beta maps from the LSA model based on condition\nCollect the :term:`Parameter Estimate` maps\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"lsa_beta_maps = {cond: [] for cond in events_df[\"trial_type\"].unique()}\ntrialwise_conditions = lsa_events_df[\"trial_type\"].unique()\nfor condition in trialwise_conditions:\n beta_map = lsa_glm.compute_contrast(condition, output_type=\"effect_size\")\n # Drop the trial number from the condition name to get the original name\n condition_name = condition.split(\"__\")[0]\n lsa_beta_maps[condition_name].append(beta_map)\n\n# We can concatenate the lists of 3D maps into a single 4D beta series for\n# each condition, if we want\nlsa_beta_maps = {\n name: image.concat_imgs(maps) for name, maps in lsa_beta_maps.items()\n}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define the LSS models\nWe will now create a separate LSS model for each trial of interest.\nThe transformation is much like the LSA approach, except that we only\nrelabel *one* trial in the DataFrame.\nWe loop through the trials, create a version of the DataFrame where the\ntargeted trial has a unique trial type, fit the model to that DataFrame,\nand finally collect the targeted trial's beta map for the beta series.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def lss_transformer(df, row_number):\n \"\"\"Label one trial for one LSS model.\n\n Parameters\n ----------\n df : pandas.DataFrame\n BIDS-compliant events file information.\n row_number : int\n Row number in the DataFrame.\n This indexes the trial that will be isolated.\n\n Returns\n -------\n df : pandas.DataFrame\n Update events information, with the select trial's trial type isolated.\n trial_name : str\n Name of the isolated trial's trial type.\n \"\"\"\n df = df.copy()\n\n # Determine which number trial it is *within the condition*\n trial_condition = df.loc[row_number, \"trial_type\"]\n trial_type_series = df[\"trial_type\"]\n trial_type_series = trial_type_series.loc[\n trial_type_series == trial_condition\n ]\n trial_type_list = trial_type_series.index.tolist()\n trial_number = trial_type_list.index(row_number)\n\n # We use a unique delimiter here (``__``) that shouldn't be in the\n # original condition names.\n # Technically, all you need is for the requested trial to have a unique\n # 'trial_type' *within* the dataframe, rather than across models.\n # However, we may want to have meaningful 'trial_type's (e.g., 'Left_001')\n # across models, so that you could track individual trials across models.\n trial_name = f\"{trial_condition}__{trial_number:03d}\"\n df.loc[row_number, \"trial_type\"] = trial_name\n return df, trial_name\n\n\n# Loop through the trials of interest and transform the DataFrame for LSS\nlss_beta_maps = {cond: [] for cond in events_df[\"trial_type\"].unique()}\nlss_design_matrices = []\n\nfor i_trial in range(events_df.shape[0]):\n lss_events_df, trial_condition = lss_transformer(events_df, i_trial)\n\n # Compute and collect beta maps\n lss_glm = FirstLevelModel(**glm_parameters)\n lss_glm.fit(fmri_file, lss_events_df)\n\n # We will save the design matrices across trials to show them later\n lss_design_matrices.append(lss_glm.design_matrices_[0])\n\n beta_map = lss_glm.compute_contrast(\n trial_condition,\n output_type=\"effect_size\",\n )\n\n # Drop the trial number from the condition name to get the original name\n condition_name = trial_condition.split(\"__\")[0]\n lss_beta_maps[condition_name].append(beta_map)\n\n# We can concatenate the lists of 3D maps into a single 4D beta series for\n# each condition, if we want\nlss_beta_maps = {\n name: image.concat_imgs(maps) for name, maps in lss_beta_maps.items()\n}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Show the design matrices for the first few trials\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(ncols=3, figsize=(20, 10))\nfor i_trial in range(3):\n plotting.plot_design_matrix(\n lss_design_matrices[i_trial],\n ax=axes[i_trial],\n )\n axes[i_trial].set_title(f\"Trial {i_trial + 1}\")\n\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare the three modeling approaches\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"DM_TITLES = [\"Standard GLM\", \"LSA Model\", \"LSS Model (Trial 1)\"]\nDESIGN_MATRICES = [\n standard_glm.design_matrices_[0],\n lsa_glm.design_matrices_[0],\n lss_design_matrices[0],\n]\n\nfig, axes = plt.subplots(\n ncols=3,\n figsize=(20, 10),\n gridspec_kw={\"width_ratios\": [1, 2, 1]},\n)\n\nfor i_ax, ax in enumerate(axes):\n plotting.plot_design_matrix(DESIGN_MATRICES[i_ax], ax=axes[i_ax])\n axes[i_ax].set_title(DM_TITLES[i_ax])\n\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Applications of beta series\nBeta series can be used much like resting-state data, though generally with\nvastly reduced degrees of freedom than a typical resting-state run, given\nthat the number of trials should always be less than the number of volumes\nin a functional MRI run.\n\nTwo common applications of beta series are to functional connectivity and\ndecoding analyses.\nFor an example of a beta series applied to decoding, see\n`sphx_glr_auto_examples_02_decoding_plot_haxby_glm_decoding.py`.\nHere, we show how the beta series can be applied to functional connectivity\nanalysis.\nIn the following section, we perform a quick task-based functional\nconnectivity analysis of each of the two task conditions\n('language' and 'string'), using the LSS beta series.\nThis section is based on\n`sphx_glr_auto_examples_03_connectivity_plot_seed_to_voxel_correlation.py`,\nwhich goes into more detail about seed-to-voxel functional connectivity\nanalyses.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n\nfrom nilearn.maskers import NiftiMasker, NiftiSpheresMasker\n\n# Coordinate taken from Neurosynth's 'language' meta-analysis\ncoords = [(-54, -42, 3)]\n\n# Initialize maskers for the seed and the rest of the brain\nseed_masker = NiftiSpheresMasker(\n coords,\n radius=8,\n detrend=True,\n standardize=\"zscore_sample\",\n low_pass=None,\n high_pass=None,\n t_r=None,\n memory=\"nilearn_cache\",\n memory_level=1,\n verbose=0,\n)\n\nbrain_masker = NiftiMasker(\n smoothing_fwhm=6,\n detrend=True,\n standardize=\"zscore_sample\",\n low_pass=None,\n high_pass=None,\n t_r=None,\n memory=\"nilearn_cache\",\n memory_level=1,\n verbose=0,\n)\n\n# Perform the seed-to-voxel correlation for the LSS 'language' beta series\nlang_seed_beta_series = seed_masker.fit_transform(lss_beta_maps[\"language\"])\nlang_beta_series = brain_masker.fit_transform(lss_beta_maps[\"language\"])\nlang_corrs = (\n np.dot(\n lang_beta_series.T,\n lang_seed_beta_series,\n )\n / lang_seed_beta_series.shape[0]\n)\nlanguage_connectivity_img = brain_masker.inverse_transform(lang_corrs.T)\n\n# Same but now for the LSS 'string' beta series\nstring_seed_beta_series = seed_masker.fit_transform(lss_beta_maps[\"string\"])\nstring_beta_series = brain_masker.fit_transform(lss_beta_maps[\"string\"])\nstring_corrs = (\n np.dot(\n string_beta_series.T,\n string_seed_beta_series,\n )\n / string_seed_beta_series.shape[0]\n)\nstring_connectivity_img = brain_masker.inverse_transform(string_corrs.T)\n\n# Show both correlation maps\nfig, axes = plt.subplots(figsize=(10, 8), nrows=2)\n\ndisplay = plotting.plot_stat_map(\n language_connectivity_img,\n threshold=0.5,\n vmax=1,\n cut_coords=coords[0],\n title=\"Language\",\n figure=fig,\n axes=axes[0],\n)\ndisplay.add_markers(\n marker_coords=coords,\n marker_color=\"g\",\n marker_size=300,\n)\n\ndisplay = plotting.plot_stat_map(\n string_connectivity_img,\n threshold=0.5,\n vmax=1,\n cut_coords=coords[0],\n title=\"String\",\n figure=fig,\n axes=axes[1],\n)\ndisplay.add_markers(\n marker_coords=coords,\n marker_color=\"g\",\n marker_size=300,\n)\nfig.suptitle(\"LSS Beta Series Functional Connectivity\")\n\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n\n .. footbibliography::\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.11.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}