# Behavioral and Decoding Analyses¶

This notebook contains all of the top-level code underlying the statistical analyses reported in the paper, although much of the heavy lifting occurs elsewhere, e.g. in the lyman package. All of the code to draw the data plots is also here.

The structure of the notebook should correspond pretty directly to the structure of the paper. All of the analysis is written in specific functions to try and keep the namespace clean and avoid bugs.

## Imports, definitions, and project-specific functions¶

First, import external packages that we'll use in this notebook

In :
from __future__ import division
import os.path as op
import itertools
import numpy as np
import scipy as sp
import pandas as pd
import nibabel as nib
from scipy import stats
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import r2_score


Now my personal libraries

In :
import lyman
from lyman import mvpa, evoked
import seaborn as sns
import moss


Set up display options for figures and tables

In :
%matplotlib inline
sns.set(context="paper", style="ticks", font="Arial")
mpl.rcParams.update({"xtick.major.width": 1, "ytick.major.width": 1, "savefig.dpi": 150})
pd.set_option('display.precision', 3)


Activate 3D plotting for matplotlib

In :
from mpl_toolkits import mplot3d
from mpl_toolkits.mplot3d import Axes3D


Set a random seed so code that uses randomness is reproducible. This mostly affects the bootstrapping that underlies error bars in the plots -- the permutation code for classifier analysis is seeded separately.

In :
np.random.seed(sum(map(ord, reversed("DKsort"))))


Load the IPython magic extension to use R within this notebook and import external R packages we'll use.

In :
%load_ext rmagic
%R library(lme4)
%R library(multcomp)

Loading required package: Matrix

Attaching package: ‘lme4’

The following objects are masked from ‘package:stats’:

AIC, BIC


Loading required package: mvtnorm


Connect to an IPython cluster (this must be started separately) and get two direct view references -- one to all engines (we expect to run this on my 8-core Mac Pro) and one to only 4 engines. We use the dv4 cluster view with the extraction functions so we can extract data in parallel without killing things by running too many parallel i/o processes.

In :
from IPython.parallel import Client, TimeoutError
try:
dv = Client()[:]
dv4 = Client()[:4]
except (IOError, TimeoutError):
dv = None
dv4 = None

/Users/mwaskom/anaconda/envs/dksort/lib/python2.7/site-packages/IPython/parallel/client/client.py:444: RuntimeWarning:
Controller appears to be listening on localhost, but not on this machine.
If this is true, you should specify Client(...,sshserver='[email protected]')
or instruct your controller to listen on an external IP.
RuntimeWarning)


Pandas 0.12 prints a very annoying deprecation warning every time an HTML representation is created for a DataFrame. We're going to silence those here:

In :
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)


Control saving figures globally with this function so it is easy to switch between file types and possibly other options.

In :
def save_figure(fig, figname):
fig.savefig("figures/%s.pdf" % figname, dpi=300)
fig.savefig("figures/%s.tiff" % figname, dpi=300)


Get the list of subject IDs for all the analyses.

In :
subjects = pd.Series(lyman.determine_subjects(), name="subj")


Get root paths for lyman outputs

In :
project = lyman.gather_project_info()
anal_dir = project["analysis_dir"]
data_dir = project["data_dir"]


Set some colors globally to use in all figures.

In :
rule_colors = dict(dim=["#EE6F25", "#4BC75B", "#9370DB"], dec=["#6495EB", "#FF6347"])
roi_colors = dict(IFS="#CC3333", aMFG="#3380CC", pMFG="#33CC80", pSFS="#CC8732", IFG="#6464D8",
aIns="#F08080", FPC="#CCCC33", IPS="#2A4380", OTC="#24913C")
roi_colors.update({"aIFS": "#863D3D", "pIFS": "#7A5252",
"lh-IFS": "#D19494", "rh-IFS": "#C2A3A3"})


Define an R function to facilitate reporting of likelihood ratio test of nested LMEMs.

In :
%%R
lr_test = function(m1, m2, name){
out = anova(m1, m2)
chi2 = out$Chisq dof = out$"Chi Df"
p = out$"Pr(>Chisq)" test_str = "Likelihood ratio test for %s:\n Chisq(%d) = %.2f; p = %.3g" writeLines(sprintf(test_str, name, dof, chi2, p)) }  ### Behavioral/design information¶ Read each subject's master behavioral file and build two group DataFrames -- one with all trials (behav_full) and one with only trials considered in the fMRI analyses (behav_df). The latter exludes incorrect trials as well as trials identified as artifacts in the fMRI preprocessing. In : def dksort_behav_join(): behav_temp = op.join(data_dir, "%s/behav/behav_data.csv") behav_data = [] for subj in subjects: df = pd.read_csv(behav_temp % subj, index_col="trial") df["subj"] = np.repeat(subj, len(df)) behav_data.append(df) behav_full = pd.concat(behav_data) behav_df = behav_full[behav_full["clean"] & behav_full["correct"]] return behav_full, behav_df behav_full, behav_df = dksort_behav_join() %Rpush behav_df %Rpush behav_full  ### fMRI Parameters and Functions¶ Define some parameters shared across fMRI analyses below. In : pfc_rois = ["IFS", "aMFG", "pMFG", "FPC", "IFG", "aIns", "pSFS"] net_rois = ["IFS", "IPS", "OTC"] all_rois = pd.Series(pfc_rois + net_rois[1:], name="roi")  In : other_rois = ["lh-IFS", "rh-IFS", "aIFS", "pIFS", "AllROIs"] other_masks = ["lh.yeo17_ifs", "rh.yeo17_ifs", "yeo17_aifs", "yeo17_pifs", "dksort_all_pfc"]  In : frames = np.arange(-1, 5) timepoints = frames * 2 + 1 up_timepoints = (np.linspace(0, 12, 25) - 1)[:-1] model = LogisticRegression() n_shuffle = 1000 peak = slice(2, 4) shuffle_seed = sum(map(ord, "DKsort"))  Define several functions to perform the same sequence of analyses on different sets of ROIs/events. In : def dksort_decode(rois, rule): """Do the all the main decoding steps across sets of ROIs. Parameters ---------- rois: list of strings list of ROI names that can be easily mapped to a mask name rule: string name of rule type corresponding to design file in data dir Returns ------- dictionary with the following entries rois: list of roi names accs: DataFrame with decoding accuracy. Index is hierarchical with (ROI, subject), and columns are timepoint. chance: DataFrame in same shape as accs with the mean value from the shuffled null distribution for each test. peak: DataFrame with decoding accuracy from data averaged over 3s and 5s. Index is subject id and column is ROI. null: DataFrame with maximum shuffled accuracy across timepoints. Index is hierarchical with (subj, iteration) and columns are ROIS. ttest: DataFrame with ROIs in the index and (t, p, max_tp) for the group average test against empirical chance in the columns. t is the t statistic for the peak accuracy value, and p is the corresponding p value corrected for multiple comparisons across region and timepoint. max_tp is in seconds relative to stimulus onset. signif: DataFrame with the null distribution percentiles corresponding to the best observed accuracy for each subject/ROI. Index is hierarchical with (correction, subject) where correction can be time or omni for the space of tests the percentile is corrected against. columns are ROIs. """ # Set up the DataFrames to hold the persisent outputs roi_index = moss.product_index([rois, subjects], ["ROI", "subj"]) columns = pd.Series(timepoints, name="timepoints") accs = pd.DataFrame(index=roi_index, columns=columns, dtype=float) null_index = moss.product_index([subjects, np.arange(n_shuffle)], ["subj", "iter"]) null = pd.DataFrame(index=null_index, columns=pd.Series(rois, name="ROI"), dtype=float) peak_df = pd.DataFrame(index=subjects, columns=pd.Series(rois, name="ROI")) chance = pd.DataFrame(index=roi_index, columns=timepoints, dtype=float) # For each ROI load the data, decode, and simulate the null distribution for roi in rois: mask = "yeo17_" + roi.lower() # Load the dataset and do the basic time-resolved decoding ds = mvpa.extract_group(rule, roi, mask, frames, confounds="rt", dv=dv4) roi_accs = mvpa.decode_group(ds, model, dv=dv) accs.loc[roi, :] = roi_accs # Now do the shuffling and save the partially transformed null distribution roi_null = mvpa.classifier_permutations(ds, model, n_iter=n_shuffle, random_seed=shuffle_seed, dv=dv) chance.loc[roi, :] = roi_null.mean(axis=1) null[roi] = roi_null.max(axis=-1).ravel() # Finally re-load the data averaging over the peak timepoints and decode peak_ds = mvpa.extract_group(rule, roi, mask, frames, peak, "rt", dv=dv4) peak_df[roi] = mvpa.decode_group(peak_ds, model, dv=dv) # Do the group t tests wide_accs = accs.unstack(level="ROI") wide_chance = chance.unstack(level="ROI") mus = wide_accs.mean(axis=0) sds = wide_accs.std(axis=0) ts, ps = moss.randomize_onesample(wide_accs, h_0=wide_chance) # Build the t test output t_df = pd.DataFrame(dict(t=ts, p=ps, mu=mus, sd=sds), index=wide_accs.columns, columns=["mu", "sd", "t", "p"]) ttest = t_df.groupby(level="ROI").apply(lambda x: x.loc[x.mu.idxmax()]) ttest["tp"] = t_df.groupby(level="ROI").mu.apply(lambda x: x.idxmax()) # From the null distribution, find the percentile corresponding to the observed score signif_index = moss.product_index([["time", "omni"], subjects], ["correction", "subj"]) signif = pd.DataFrame(index=signif_index, columns=rois, dtype=float) for roi, roi_accs in accs.max(axis=1).groupby(level="ROI"): for subj, score in roi_accs.groupby(level="subj"): signif.loc[("time", subj), roi] = stats.percentileofscore(null[roi], score) signif.loc[("omni", subj), roi] = stats.percentileofscore(null.max(axis=1), score) return dict(rois=rois, accs=accs, chance=chance, peak=peak_df, null=null, ttest=ttest, signif=signif)  In : def dksort_timecourse_figure(ax, data, ytick_args, err_style="ci_band", legend=True, err_kws=None, **kwargs): """Plot time-resolved decoding accuracies for multiple ROIs""" # Represent chance empirically based on the null distribution chance = np.array(data["chance"].mean(axis=0)) ax.plot(timepoints, chance, "k--") # Plot the stimulus onset ax.plot([0, 0], ytick_args[:2], ls=":", c="k") # Draw the accuracy timecourse for each ROI colors = [roi_colors[roi] for roi in data["rois"]] interpolate = not err_style == "ci_bars" accs = pd.melt(data["accs"].reset_index(), ["subj", "ROI"], value_name="acc") sns.tsplot(accs, time="timepoints", unit="subj", condition="ROI", value="acc", color=colors, err_style=err_style, interpolate=interpolate, legend=legend, err_kws=err_kws, ax=ax, **kwargs) # Adjust the axis scales and tick labels ylim = ytick_args[:2] yticks = np.linspace(*ytick_args) ax.set_xlim(timepoints.min(), timepoints.max()) ax.set_ylim(ylim) ax.set_yticks(yticks) # Label the axes ax.set_xlabel("Time relative to stimulus onset (s)") ax.set_ylabel("Cross-validated decoding accuracy") # Make a legend with the ROIs ax.legend(frameon=False, loc="upper right")  In : def dksort_point_figure(ax, data, chance, ytick_args, yaxis_label): """Plot a point estimate and CI for the 3-5s accuracy by ROI.""" # Plot the chance line ax.axhline(chance, color="k", ls="--") # Plot each accuracy and confidence interval for i, roi in enumerate(data.columns): color = roi_colors[roi] accs = data[roi].values ci = moss.ci(moss.bootstrap(accs), 68) ax.plot(i, accs.mean(), "o", color=color, ms=5, mew=0, mec=color) ax.plot([i, i], ci, color=color, lw=2, solid_capstyle="round") # Set the axis limits ax.set_xlim(-.5, data.shape - .5) ax.set_ylim(*ytick_args[:2]) yticks = np.linspace(*ytick_args) ax.set_yticks(yticks) ax.set_yticklabels(["%.2f" % t for t in yticks]) ax.xaxis.grid(False) # Set the axis labels ax.set_xticks(np.arange(len(data.columns))) ax.set_xticklabels(data.columns, rotation=60) if yaxis_label: ax.set_ylabel("Cross-validated decoding accuracy") else: ax.set_yticklabels([])  In : def dksort_significance_figure(ax, data): """Plot the percentile of the null for each ROI by subject.""" width = 2 / len(data) xbase = np.arange(0, 1, width) rois = data.columns # Plot percentile across ROI/time and time correction for kind, kind_df in data.groupby(level="correction"): if kind == "time": continue for i, roi in enumerate(rois): height = kind_df[roi] color = sns.desaturate(roi_colors[roi], .75) edgecolor = sns.desaturate(color, .5) alpha = .25 if kind == "time" else .7 ax.bar(xbase + i, np.sort(height), width=width, color=color, edgecolor=edgecolor, lw=0.1, alpha=alpha) if i: ax.axvline(i, lw=.5, color="#555555") # Set the other plot aesthetics ax.set_ylim(0, 100) ax.axhline(95, c="#444444", ls=":") ax.set_xticks(np.arange(data.shape) + .5) ax.set_xticklabels(rois, rotation=60) ax.xaxis.grid(False) ax.set_ylabel("Percentile in null distribution")  ## Behavioral results¶ ### Speed and accuracy summary¶ Summarize overall accuracy and reaction time across subjects. Take medians for RT within subject, but report the group mean over these medians. In : def dksort_behav_summary(): accs = behav_full.groupby("subj").correct.mean() print "Accuracy across subjects: Mean %.2g; Std %.2g; Min: %.2g" % (accs.mean(), accs.std(), accs.min()) rts = behav_df.groupby("subj").rt.median() * 1000 print "Median RT across subjects (ms) - Mean %.3g; Std %.3g" % (rts.mean(), rts.std()) iqrs = behav_df.groupby("subj").rt.apply(moss.iqr) * 1000 print "RT IQRs across subjects (ms) - Mean: %.3g; Std %.3g" % (iqrs.mean(), iqrs.std()) dksort_behav_summary()  Accuracy across subjects: Mean 0.95; Std 0.033; Min: 0.87 Median RT across subjects (ms) - Mean 826; Std 209 RT IQRs across subjects (ms) - Mean: 307; Std 101  ### Speed as a function of rule¶ Next, test whether the task rules influence the speed of responses. We'll perform these analyses in R using linear mixed effects models (LMEMs), which account for both within- and between-subjects variance. To test the significance of the models, we'll use likelihood ratio tests (Barr et al. 2013). Also following this paper, we will use random effects structures that model random slopes for all main effect terms considered in the fixed effects structure. First, we fit both an interactive model and additive model (betwen rulesets) predicting RT and use a likelihood ratio test to see whether the interaction is significant. In : %%R m.rt.int = lmer(rt ~ dim_rule * dec_rule + (dim_rule + dec_rule | subj), behav_df, REML=F) m.rt.add = lmer(rt ~ dim_rule + dec_rule + (dim_rule + dec_rule | subj), behav_df, REML=F) lr_test(m.rt.int, m.rt.add, "interaction")  Likelihood ratio test for interaction: Chisq(2) = 0.57; p = 0.752  Print the additive model coefficients and standard errors In : %R print(m.rt.add, corr=FALSE)  Linear mixed model fit by maximum likelihood Formula: rt ~ dim_rule + dec_rule + (dim_rule + dec_rule | subj) Data: behav_df AIC BIC logLik deviance REMLdev 925.4 1020 -447.7 895.4 920.1 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 0.04205506 0.205073 dim_rulepattern 0.00140533 0.037488 -0.350 dim_ruleshape 0.00094827 0.030794 -0.130 0.974 dec_rulesame 0.00291850 0.054023 -0.022 0.174 0.179 Residual 0.07148692 0.267370 Number of obs: 3961, groups: subj, 15 Fixed effects: Estimate Std. Error t value (Intercept) 0.923318 0.053624 17.218 dim_rulepattern -0.009415 0.014222 -0.662 dim_ruleshape -0.016024 0.013077 -1.225 dec_rulesame -0.083386 0.016337 -5.104  Report the omnibus test for the additive model and likelihood ratio tests for each main effect. In : %%R print(anova(m.rt.add)) m.rt.dim = lmer(rt ~ dec_rule + (dim_rule + dec_rule | subj), behav_df, REML=F) m.rt.dec = lmer(rt ~ dim_rule + (dim_rule + dec_rule | subj), behav_df, REML=F) lr_test(m.rt.add, m.rt.dim, "dimension rule") lr_test(m.rt.add, m.rt.dec, "decision rule")  Analysis of Variance Table Df Sum Sq Mean Sq F value dim_rule 2 0.05954 0.02977 0.4164 dec_rule 1 1.86241 1.86241 26.0524 Likelihood ratio test for dimension rule: Chisq(2) = 1.52; p = 0.467 Likelihood ratio test for decision rule: Chisq(1) = 15.09; p = 0.000103  It's possible that the dimension rules may influence RTs differently for each participant and wash out in the group analysis. To examine this possiblitly, perform a one-way ANOVA across each ruleset within each participant and count the number of subjects with a trending or significant effect. In : def dksort_rt_anova(): subj_grouped = behav_df.groupby("subj") dim = subj_grouped.apply(moss.df_oneway, "dim_rule", "rt", False) dec = subj_grouped.apply(moss.df_oneway, "dec_rule", "rt", False) f_scores = pd.concat([dim, dec], keys=["dimension", "decision"], names=["rules", "subj"]) print "Median F score: %.2f" % f_scores.loc["dimension"]["F"].median() f_scores["p < 0.05"] = f_scores.p < 0.05 f_scores["p < 0.1"] = f_scores.p < 0.1 print f_scores.groupby(level="rules")[["p < 0.05", "p < 0.1"]].sum() return subj_grouped, f_scores subj_grouped, f_scores = dksort_rt_anova()  Median F score: 0.67 p < 0.05 p < 0.1 rules dimension 2 2 decision 8 10  Also, to get a measurement of the slowing associated with the 'different' rule, take the difference in medians across the rules within subjects, then report the mean and bootstrapped confidence interval for the mean across subjects. In : def dksort_rule_rt_effect(): pivot = pd.pivot_table(behav_df, "rt", "subj", "dec_rule", np.median) cost = np.diff(pivot, axis=1) * 1000 boots = moss.bootstrap(cost) ci_low, ci_high = moss.ci(boots, 95) print "'Different' rule effect (ms): Mean %.3g, CI: [%.3g, %.3g]" % (boots.mean(), ci_low, ci_high) dksort_rule_rt_effect()  'Different' rule effect (ms): Mean -94.4, CI: [-123, -66.2]  ### Speed as a function of attended feature matching and decision rule¶ Next we'll ask whether RT is modulated by whether the attended features match In : %%R m.match = lmer(rt ~ attend_match + (attend_match | subj), behav_df, REML=FALSE) m.match.drop = lmer(rt ~ (attend_match | subj), behav_df, REML=FALSE) print(m.match, corr=FALSE) lr_test(m.match, m.match.drop, "attended matching")  Linear mixed model fit by maximum likelihood Formula: rt ~ attend_match + (attend_match | subj) Data: behav_df AIC BIC logLik deviance REMLdev 1028 1066 -508.2 1016 1028 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 0.04046794 0.201166 attend_matchTRUE 0.00014153 0.011897 0.208 Residual 0.07422838 0.272449 Number of obs: 3961, groups: subj, 15 Fixed effects: Estimate Std. Error t value (Intercept) 0.866011 0.052302 16.558 attend_matchTRUE 0.013382 0.009189 1.456 Likelihood ratio test for attended matching: Chisq(1) = 1.99; p = 0.159  We'll also ask whether this interacts with the decision rule. In : %%R m.match.int = lmer(rt ~ attend_match * dec_rule + (attend_match + dec_rule | subj), behav_df, REML=FALSE) m.match.add = lmer(rt ~ attend_match + dec_rule + (attend_match + dec_rule | subj), behav_df, REML=FALSE) print(m.match.int, corr=FALSE) lr_test(m.match.int, m.match.add, "interaction")  Linear mixed model fit by maximum likelihood Formula: rt ~ attend_match * dec_rule + (attend_match + dec_rule | subj) Data: behav_df AIC BIC logLik deviance REMLdev 842.3 911.4 -410.1 820.3 844.7 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 0.03968419 0.199209 attend_matchTRUE 0.00017723 0.013313 0.322 dec_rulesame 0.00299747 0.054749 0.010 -0.782 Residual 0.07028425 0.265112 Number of obs: 3961, groups: subj, 15 Fixed effects: Estimate Std. Error t value (Intercept) 0.870299 0.052129 16.695 attend_matchTRUE 0.088947 0.012443 7.148 dec_rulesame -0.007847 0.018508 -0.424 attend_matchTRUE:dec_rulesame -0.150621 0.016855 -8.936 Likelihood ratio test for interaction: Chisq(1) = 79.04; p = 6.08e-19  ### Response accuracy as a function of rule¶ Now perform a similar sequence of analyses on response accuracy, here using mixed effects logistic regression models. In : %%R m.acc.int = lmer(correct ~ dim_rule * dec_rule + (dim_rule + dec_rule | subj), behav_full, family=binomial) m.acc.add = lmer(correct ~ dim_rule + dec_rule + (dim_rule + dec_rule | subj), behav_full, family=binomial) lr_test(m.acc.int, m.acc.add, "interaction")  Likelihood ratio test for interaction: Chisq(2) = 0.15; p = 0.927  In : %R print(m.acc.add, corr=FALSE)  Generalized linear mixed model fit by the Laplace approximation Formula: correct ~ dim_rule + dec_rule + (dim_rule + dec_rule | subj) Data: behav_full AIC BIC logLik deviance 1696 1785 -833.8 1668 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 0.46217 0.67983 dim_rulepattern 0.01600 0.12649 -1.000 dim_ruleshape 0.22851 0.47802 -0.250 0.250 dec_rulesame 0.15589 0.39483 -0.285 0.285 -0.347 Number of obs: 4320, groups: subj, 15 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.0981 0.2304 13.449 <2e-16 *** dim_rulepattern -0.3197 0.1775 -1.801 0.0716 . dim_ruleshape -0.1417 0.2229 -0.636 0.5249 dec_rulesame 0.4053 0.1805 2.246 0.0247 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  In : %%R m.acc.dim = lmer(correct ~ dec_rule + (dim_rule + dec_rule | subj), behav_full, family=binomial) m.acc.dec = lmer(correct ~ dim_rule + (dim_rule + dec_rule | subj), behav_full, family=binomial) lr_test(m.acc.add, m.acc.dim, "dimension rules") lr_test(m.acc.add, m.acc.dec, "decision rules")  Likelihood ratio test for dimension rules: Chisq(2) = 2.64; p = 0.267 Likelihood ratio test for decision rules: Chisq(1) = 3.85; p = 0.0498  In : def dksort_rule_acc_effect(): pivot = pd.pivot_table(behav_full, "correct", "subj", "dec_rule") cost = np.diff(pivot, axis=1) * 100 boots = moss.bootstrap(cost) low, high = moss.ci(boots, 95) args = cost.mean(), low, high print "'Different' rule response accuracy cost: %.2g%%; CI: [%.2g%%, %.2g%%]" % args dksort_rule_acc_effect()  'Different' rule response accuracy cost: 1.9%; CI: [0.28%, 3.7%]  ### Rule Switching Analyses¶ Now consider how shifting back and forth between different rulesets influences response speed. First, analyze it as a function of how many trials have elapsed since the rule switch. We'll fit this in two separate models as the lag is frequently identical for each ruleset. In : %%R m.lag.dim = lmer(rt ~ dim_shift_lag + (dim_shift_lag | subj), behav_df, REML=FALSE) m.lag.dec = lmer(rt ~ dec_shift_lag + (dec_shift_lag | subj), behav_df, REML=FALSE)  In : %R print(m.lag.dim, corr=FALSE)  Linear mixed model fit by maximum likelihood Formula: rt ~ dim_shift_lag + (dim_shift_lag | subj) Data: behav_df AIC BIC logLik deviance REMLdev 1031 1068 -509.3 1019 1033 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 4.0091e-02 0.2002269 dim_shift_lag 1.0637e-06 0.0010314 1.000 Residual 7.4303e-02 0.2725866 Number of obs: 3961, groups: subj, 15 Fixed effects: Estimate Std. Error t value (Intercept) 0.8714653 0.0520818 16.733 dim_shift_lag 0.0005944 0.0022995 0.258  In : %R print(m.lag.dec, corr=FALSE)  Linear mixed model fit by maximum likelihood Formula: rt ~ dec_shift_lag + (dec_shift_lag | subj) Data: behav_df AIC BIC logLik deviance REMLdev 1022 1059 -504.9 1010 1024 Random effects: Groups Name Variance Std.Dev. Corr subj (Intercept) 4.4698e-02 0.2114198 dec_shift_lag 1.5458e-05 0.0039316 -0.830 Residual 7.4109e-02 0.2722290 Number of obs: 3961, groups: subj, 15 Fixed effects: Estimate Std. Error t value (Intercept) 0.885098 0.054989 16.096 dec_shift_lag -0.004589 0.002155 -2.129  In : %%R m.lag.nodim = lmer(rt ~ 1 + (dim_shift_lag | subj), behav_df, REML=FALSE) m.lag.nodec = lmer(rt ~ 1 + (dec_shift_lag | subj), behav_df, REML=FALSE) lr_test(m.lag.dim, m.lag.nodim, "dimension rules") lr_test(m.lag.dec, m.lag.nodec, "decision rules")  Likelihood ratio test for dimension rules: Chisq(1) = 0.07; p = 0.797 Likelihood ratio test for decision rules: Chisq(1) = 3.88; p = 0.0487  Next, directly test switch costs by analyzing RT on the first trial of each block as a function of rule switches. We can do this in one model. In : %%R shift_subset = behav_df$block_pos == 0
m.shift.add = lmer(rt ~ dim_shift + dec_shift + (dim_shift + dec_shift | subj),
behav_df, subset=shift_subset, REML=FALSE)
m.shift.int = lmer(rt ~ dim_shift * dec_shift + (dim_shift + dec_shift | subj), behav_df,
subset=shift_subset, REML=FALSE)

Likelihood ratio test for interaction:
Chisq(1) = 0.61; p = 0.436

In :
%R print(m.shift.add, corr=FALSE)

Linear mixed model fit by maximum likelihood
Formula: rt ~ dim_shift + dec_shift + (dim_shift + dec_shift | subj)
Data: behav_df
Subset: shift_subset
AIC   BIC logLik deviance REMLdev
397 448.6 -188.5      377   393.7
Random effects:
Groups   Name          Variance   Std.Dev. Corr
subj     (Intercept)   0.04266668 0.206559
dim_shiftTRUE 0.00013579 0.011653 -1.000
dec_shiftTRUE 0.00056459 0.023761  1.000 -1.000
Residual               0.07488755 0.273656
Number of obs: 1290, groups: subj, 15

Fixed effects:
Estimate Std. Error t value
(Intercept)    0.86474    0.05606  15.426
dim_shiftTRUE  0.01119    0.01728   0.648
dec_shiftTRUE  0.01977    0.01668   1.185

In :
%%R
m.shift.nodim = lmer(rt ~ dec_shift + (dim_shift + dec_shift | subj),
behav_df, subset=shift_subset, REML=FALSE)
m.shift.nodec = lmer(rt ~ dim_shift + (dim_shift + dec_shift | subj),
behav_df, subset=shift_subset, REML=FALSE)

Likelihood ratio test for dimension rule:
Chisq(1) = 0.42; p = 0.517
Likelihood ratio test for decision rule:
Chisq(1) = 1.39; p = 0.239


### Behavioral analysis figures¶

Here is the main behavioral analysis figure. There is a lot of information packed into this plot, so it is complicated to draw.

In :
def dksort_figure_2():
# Set up variables shared across subplots
dfs = [behav_df, behav_full]
measures = ["rt", "correct"]
agg_funcs = [np.median, np.mean]
ci = 68
xlabels = ["Dimension rules", "Attended features", "Trials since rule switch"]
ylabels = ["Reaction time (s)", "Response accuracy"]
xticks = [range(3), range(2), range(1, 7)]
xticklabels = [["Shape", "Color", "Pattern"], ["Mismatch", "Match"], range(1, 7)]
xlims = [[-.5, 2.5], [-.5, 1.5], [.25, 6.75]]
ylims = [[.67, 1.02], [.8, 1.02]]
yticks = [[.7, 1, 4], [.8, 1, 3]]
ms, mew, lw = 3.5, 1.2, 1.2
err_kws = dict(linewidth=lw, mew=mew)
lag_colors = dict(Dimension="black", Decision="lightslategray")
text_offset = (0.01, -0.01)
text_size = 11

# Draw each plot
f, axes = plt.subplots(2, 3, figsize=(4.48, 3))
for i in range(2):

# First analyze RT/accuracy analysis sorted by rule types
pivots = []
for subj, df_subj in dfs[i].groupby("subj"):
pivot = pd.pivot_table(df_subj, measures[i], "dim_rule", "dec_rule", agg_funcs[i])
pivots.append(np.array(pivot))
pivots = np.array(pivots, float)
means = pivots.mean(axis=0)
cis = moss.ci(moss.bootstrap(pivots, axis=0), ci, axis=0)

# Plot the above results in plot columns
ax = axes[i, 0]
for dim in range(3):
color = rule_colors["dim"][dim]
for dec in range(2):
x = dim + (dec - .5) / 3.5
y = means[dim, dec]
ci_x = cis[:, dim, dec]
mfc = ax.get_axis_bgcolor() if dec else color

# Plotting calls
ax.plot([x, x], ci_x, lw=lw, color=color)
ax.plot(x, y, "o", ms=ms, mew=mew, mfc=mfc, mec=color)

# Draw two plots out of the axis range to support the legend
for j, rule in enumerate(["Same", "Different"]):

# Plotting calls
ax.plot(-1, -1, "o", color="#444444", mec="#444444", ms=ms,
mfc=["none", "#444444"][j], mew=mew, label="'%s' rule" % rule)

# Next plot the decision rule by feature matching analysis
pivots = []
for subj, df_subj in dfs[i].groupby("subj"):
pivot = pd.pivot_table(df_subj, measures[i], "attend_match", "dec_rule", agg_funcs[i])
pivots.append(np.array(pivot))
pivots = np.array(pivots, float)
means = pivots.mean(axis=0)
cis = moss.ci(moss.bootstrap(pivots, axis=0), ci, axis=0)

ax  = axes[i, 1]
for k, rule in enumerate(["Same", "Different"]):
color = rule_colors["dec"][k]
x = np.array([0, 1])

# Plotting calls
sns.tsplot(pivots[..., 1 - k], time=x, err_style="ci_bars", ci=ci,
color=color, marker="o", mec=color, label = "'%s' rule" % rule,
ms=ms, mfc=color, mew=mew, lw=lw, err_kws=err_kws, ax=ax)

# Now do the analysis as a function of rule switching
for j, ruleset in enumerate(["Dimension", "Decision"]):
color = lag_colors[ruleset]
pivot_cols = ruleset[:3].lower() + "_shift_lag"
lag_pivot = pd.pivot_table(dfs[i], measures[i], "subj", pivot_cols, agg_funcs[i])
ts_kws =  dict(err_style="ci_bars", color=color, ci=ci,
lw=lw, err_kws=err_kws, ax=axes[i, 2], marker="o", ms=3)

# Plotting calls
sns.tsplot(lag_pivot.loc[:, :2].values, time=[1, 2, 3], **ts_kws)
sns.tsplot(lag_pivot.loc[:, (2, 3)].values, time=[3, 4], linestyle=":", **ts_kws)
sns.tsplot(lag_pivot.loc[:, 3:5].values, time=[4, 5, 6], label=ruleset + " rules", **ts_kws)

# Finally, iterate through the plots and set more supporting variables
for j in range(3):
ax = axes[i, j]
yticks_ = np.linspace(*yticks[i])
ax.set_xticks(xticks[j])
ax.set_yticks(yticks_)
ax.set_xlim(*xlims[j])
ax.set_ylim(ylims[i])
if i:
ax.set_xlabel(xlabels[j])
ax.set_xticklabels(xticklabels[j])
else:
ax.set_xticklabels([])
if not j:
ax.set_ylabel(ylabels[i])
ax.set_yticklabels(yticks_)
else:
ax.set_yticklabels([])
if i:
ax.legend(loc="lower left")

f.subplots_adjust(wspace=.08, hspace=.08, left=.1, bottom=.13, top=.97, right=.97)

# Label the facets
for ax, s in zip(f.axes, "ACEBDF"):
(x, _), (_, y) = ax.bbox.transformed(f.transFigure.inverted()).get_points()
x += text_offset
y += text_offset
f.text(x, y, s, size=text_size, ha="left", va="top")

sns.despine()
save_figure(f, "figure_2")

dksort_figure_2() **Main behavioral results.** Mean within-subject median reaction time (top row) and mean within-subject proportion correct responses (bottom row) sorted by the two rule sets (A and B), by the decision rule and whether the attended features matched or differed (C and D), and by the number of since the rule switch plotted separately with respect to the dimension and decision rules (E and F). The timecourses in E and F are dashed between trials 3 and 4 to indicate passage through a mini-block boundary. Error bars on all facets indicate bootstrapped standard errors of the aggregate values across subjects. Reaction times are plotted only for trials included in the imaging analyses.

Now the supplemental behavioral figure (not used in the JNeuro submission, but keeping this around).

In :
def dksort_supplemental_figure_1():
f, axes = plt.subplots(2, 1, figsize=(6.85, 6.3))

positions = dict(dim=[.25, .5, .75], dec=[.33, .66])
widths = dict(dim=.18, dec=.25)
rules = dict(dim=["shape", "color", "pattern"], dec=["same", "different"])

text_offset = (.01, -.015)
text_size = 12

subj_medians = subj_grouped.rt.median()
subj_sorted = subjects[np.argsort(subj_medians)]
for ruleset, ax in zip(["dimension", "decision"], axes):
rs = ruleset[:3]
for i, subj in enumerate(subj_sorted):
df_i = behav_df[behav_df["subj"] == subj]
for j, rule in enumerate(rules[rs]):
df_rule = df_i[df_i[rs + "_rule"] == rule]
pos = np.array([i + positions[rs][j]])
rule_rt = np.array(df_rule["rt"])
color = sns.desaturate(rule_colors[rs][j], .7)
sns.boxplot([rule_rt], color=color, positions=pos, fliersize=2,
whis=1.75, linewidth=1, widths=widths[rs], ax=ax)
ax.set_xlim(0, 15)
ax.set_xticks(np.arange(15) + 0.5)
rule_fs = f_scores["F"][ruleset].reindex(subj_sorted)
ax.set_xticklabels(["%.2g" % _f for _f in rule_fs])
ax.set_ylabel("Reaction time (s)")
ax.set_xlabel("One-way F score over rules")
indiv_rules = rules[ruleset[:3]]
n_rules = len(indiv_rules)
ax.legend(ax.artists[:n_rules], indiv_rules, loc=4,
ncol=n_rules)
sns.despine()
f.tight_layout()

f.text(.02, .96, "A", size=text_size)
f.text(.02, .47, "B", size=text_size)

save_figure(f, "supplemental_figure_1")

dksort_supplemental_figure_1() ## Descriptive statistics about ROIs and artifacts¶

### ROI sizes¶

For each subject and ROI, load the mask file and count the number of voxels included in the mask.

In :
def dksort_roi_sizes():
roi_sizes = pd.DataFrame(columns=all_rois, index=subjects, dtype=float)
for roi in roi_sizes:
for subj in subjects:
roi_sizes.loc[subj, roi] = vox_count
return roi_sizes.describe().loc[["mean", "std"]].astype(int)
dksort_roi_sizes()

Out:
roi IFS aMFG pMFG FPC IFG aIns pSFS IPS OTC
mean 626 573 468 595 534 356 277 409 1865
std 83 77 68 70 63 37 61 52 218

### Mean signal¶

Now look at the mean signal across each ROI in the original (preprocessed) timecourse.

In :
def dksort_roi_signal():
roi_signal = pd.DataFrame(columns=all_rois, index=subjects, dtype=float)
for roi in all_rois:
signal = evoked.extract_group("yeo17_" + roi.lower(), dv=dv4)
signal = [np.concatenate(d["data"]).mean() for d in signal]
roi_signal[roi] = signal
return (roi_signal / 100).describe().loc[["mean", "std"]]
dksort_roi_signal()

Out:
roi IFS aMFG pMFG FPC IFG aIns pSFS IPS OTC
mean 116.36 118.42 117.22 101.45 99.08 117.86 116.66 138.02 117.85
std 2.51 5.29 4.75 7.78 6.01 2.75 3.80 5.81 5.37

### Functional artifacts¶

In :
def dksort_artifact_counts():
artifacts = []
for subj in subjects:
subj_artifacts = 0
for run in range(1, 5):
art = pd.read_csv(op.join(anal_dir, "dksort/%s/preproc/run_%d/artifacts.csv" % (subj, run)))
art = art.max(axis=1)
subj_artifacts += art.sum()
artifacts.append(subj_artifacts)
print "Mean number of artifacts: %d" % np.mean(artifacts)
print "Percent artifact scans: %.1f%%" % (np.mean(artifacts) / 2328 * 100)
dksort_artifact_counts()

Mean number of artifacts: 23
Percent artifact scans: 1.0%


## Context decoding in lateral PFC¶

Now we turn to the central fMRI analyses: multivariate decoding of task variables from voxel populations in our cortical ROIs.

### Decoding of dimension rules¶

First try to decode the Dimension Rules -- that is, whether subjects should attend to the color, shape, or pattern of the multidimensional stimuli. We will fit separate models for each of the six prefrontal ROIs at each of six timepoints from one TR before stimulus onset (during which time there was an anticipatory cue that signlaed only an oncoming stimulus) through five TRs following stimulus onset.

We defined a function above that does the heavy lifting on the analysis. We'll run that and then unpack the results over the next few cells.

In :
pfc_dimension = dksort_decode(pfc_rois, "dimension")


We calculate a t statistic for each region and timepoint by subtracting a shuffled chance estimate from each subject's accuracy and then using a randomization test against the null hypothesis that the group mean is 0. We then report the t statiatisc for the peak accuracy across time for each region, along with a p value that is corrected for multiple comparisons across time and region.

In :
pfc_dimension["ttest"]

Out:
mu sd t p tp
ROI
FPC 0.36 0.03 3.31 7.12e-02 7
IFG 0.36 0.04 2.44 2.96e-01 3
IFS 0.46 0.06 7.95 1.00e-04 3
aIns 0.37 0.03 4.48 8.70e-03 3
aMFG 0.39 0.04 6.60 5.00e-04 5
pMFG 0.38 0.02 8.69 0.00e+00 5
pSFS 0.38 0.03 5.92 7.00e-04 5

We also use the shuffled null distribution to assess the significance of the individual subject decoding models. Report here how many subjects have at least one model over time for each region reach "significance", where signficiance is defined as p < 0.05, one-tailed, corrected for multiple comparisons either across regions and time or only across time.

In :
pfc_dimension["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())

Out:
IFS aMFG pMFG FPC IFG aIns pSFS
correction
omni 11 2 0 3 1 1 1
time 13 4 6 4 4 3 4

Perform a paired T-test on the "peak" decoding accuracy between each PFC region. Use a Bonferroni correction on the P values for the 15-way comparisons.

We don't actually use this figure, but the function presents a nice table.

In :
def dksort_pfc_paired_ttests():
pfc_paired_t = pd.DataFrame(index=pfc_rois, columns=pfc_rois, dtype=float)
pfc_paired_p = pd.DataFrame(index=pfc_rois, columns=pfc_rois, dtype=float)
for roi_i, roi_j in itertools.product(pfc_rois, pfc_rois):
peak_i = pfc_dimension["peak"][roi_i]
peak_j = pfc_dimension["peak"][roi_j]
t, p = stats.ttest_rel(peak_i, peak_j)
pfc_paired_t.loc[roi_i, roi_j] = t
pfc_paired_p.loc[roi_i, roi_j] = p
pfc_paired_p *= 15
sns.symmatplot(pfc_paired_t.values, pfc_paired_p.values,
names=pfc_rois, cmap="RdPu_r", cmap_range=(-10, 0))

dksort_pfc_paired_ttests() ### Dimension rule decoding in PFC with other ROI definitions¶

The next set of analyses support the localization claims about dimension context representation in the IFS. This is a bit of a hodepodge, but we a) break up the IFS two ways (lateralization and anterior/posterior division), b) test the IFS against an agglomerative region with all PFC ROIs together, and c) test decoding after averaging the IFS signal across the voxels.

We do all of these analyses on the "peak" data, which averages the BOLD data between the 3s and 5s timepoints before fitting the decoding models.

In :
def dksort_peak_decode(problem, rois, masks, comparison):
rois = pd.Series(rois, name="ROI")
accs = pd.DataFrame(columns=rois, index=subjects, dtype=np.float)
ttest = pd.DataFrame(columns=rois, index=["t", "p"],  dtype=np.float)
ds = mvpa.extract_group(problem, roi, mask, frames, peak, "rt", dv=dv4)
roi_accs = mvpa.decode_group(ds, model, dv=dv).squeeze()
accs[roi] = roi_accs
ttest[roi] = stats.ttest_rel(roi_accs, comparison)
return dict(accs=accs, ttest=ttest, descrip=accs.describe())

In :
ifs_subrois = dksort_peak_decode("dimension", other_rois, other_masks, pfc_dimension["peak"]["IFS"])

In :
def dksort_ifs_mean_decode():
mean_ds = mvpa.extract_group("dimension", "IFS", "yeo17_ifs", frames, peak, "rt")

for ds in mean_ds:
ds["X"] = ds["X"].mean(axis=1).reshape(-1, 1)
ds["roi_name"] = "IFS_mean"

accs = mvpa.decode_group(mean_ds, model, dv=dv).squeeze()
null = mvpa.classifier_permutations(mean_ds, model, n_shuffle,
random_seed=shuffle_seed, dv=dv).squeeze()
thresh = moss.percentiles(null, 95, axis=1)

ifs_subrois["accs"]["IFS_mean"] = accs
ifs_subrois["ttest"]["IFS_mean"] = stats.ttest_rel(accs, pfc_dimension["peak"]["IFS"])

return dict(accs=accs, null=null, thresh=thresh)

ifs_mean = dksort_ifs_mean_decode()


Look at the summary statistics for these ROIs across subjects.

In :
ifs_subrois["accs"].describe().T

Out:
count mean std min 25% 50% 75% max
ROI
lh-IFS 15 0.45 0.05 0.37 0.41 0.43 0.48 0.57
rh-IFS 15 0.44 0.05 0.38 0.39 0.45 0.48 0.52
aIFS 15 0.45 0.06 0.35 0.40 0.45 0.50 0.55
pIFS 15 0.45 0.04 0.39 0.43 0.45 0.47 0.53
AllROIs 15 0.44 0.05 0.36 0.40 0.44 0.47 0.54
IFS_mean 15 0.35 0.03 0.30 0.32 0.36 0.37 0.39
In :
ifs_subrois["accs"].apply(lambda d: pd.Series(moss.ci(moss.bootstrap(d), 95), index=[2.5, 97.5]), axis=0)

Out:
ROI lh-IFS rh-IFS aIFS pIFS AllROIs IFS_mean
2.5 0.42 0.42 0.42 0.43 0.41 0.34
97.5 0.47 0.47 0.48 0.47 0.47 0.36

Look at the paired T test for each "other" ROI against the original IFS data.

In :
ifs_subrois["ttest"].T

Out:
t p
ROI
lh-IFS -4.25 8.11e-04
rh-IFS -5.31 1.10e-04
aIFS -3.82 1.87e-03
pIFS -4.14 1.01e-03
AllROIs -4.71 3.33e-04
IFS_mean -9.24 2.46e-07

Directly test for laterality.

In :
def dksort_ifs_lateralization_test():
diff = ifs_subrois["accs"]["lh-IFS"] - ifs_subrois["accs"]["rh-IFS"]
low, high = moss.ci(moss.bootstrap(diff), 95)
t, p = stats.ttest_1samp(diff, 0)
args = (diff.mean(), low, high, t, p)
print "Test for laterality: mean difference = %.3f; 95%% CI: %.3f, %.3f, t = %.2f; p = %.3g " % args

dksort_ifs_lateralization_test()

Test for laterality: mean difference = 0.004; 95% CI: -0.023, 0.028, t = 0.29; p = 0.775


Directly test for a rostro-caudal difference in decoding performance.

In :
def dksort_ifs_rostral_test():
diff = ifs_subrois["accs"]["aIFS"] - ifs_subrois["accs"]["pIFS"]
low, high = moss.ci(moss.bootstrap(diff), 95)
t, p = stats.ttest_1samp(diff, 0)
print "Test for rosto-caudal organization: mean difference = %.3f; " % diff.mean(),
print "95%% CI: %.3f, %.3f, t = %.2f; p = %.3g " % (low, high, t, p)

dksort_ifs_rostral_test()

Test for rosto-caudal organization: mean difference = 0.002;  95% CI: -0.022, 0.024, t = 0.13; p = 0.896


Test the reduced dimensionality dataset against chance.

In :
def dksort_ifs_mean_test():
low, high = moss.ci(moss.bootstrap(ifs_mean["accs"]), 95)
chance = ifs_mean["null"].mean()
delta = ifs_mean["accs"] - chance
t, p = moss.randomize_onesample(delta)
signif = ifs_mean["accs"] > ifs_mean["thresh"]

print "IFS_Mean ROI:"
print "  mean accuracy = %.3f; 95%% CI: %.3f, %.3f;" % (ifs_mean["accs"].mean(), low, high),
print "chance = %.3f; t = %.2f; p = %.3g " % (chance, t, p)
print "  %d/15 significant IFS_Mean models" % signif.sum()

dksort_ifs_mean_test()

IFS_Mean ROI:
mean accuracy = 0.350; 95% CI: 0.336, 0.364; chance = 0.339; t = 1.52; p = 0.0729
2/15 significant IFS_Mean models


Now we have all of the information we need to make the figure about decoding dimension rules in PFC

In :
def dksort_figure_4():
f = plt.figure(figsize=(4.48, 5.5))

ax_ts = f.add_axes([.12, .53, .86, .45])
dksort_timecourse_figure(ax_ts, pfc_dimension, (.32, .48, 5))

ax_sig = f.add_axes([.12, .08, .34, .34])
dksort_significance_figure(ax_sig, pfc_dimension["signif"])

ax_pfc_peak = f.add_axes([.59, .08, .22, .34])
dksort_point_figure(ax_pfc_peak, pfc_dimension["peak"], .33, (.32, .5, 7), True)

ax_sub_peak = f.add_axes([.86, .08, .12, .34])
dksort_point_figure(ax_sub_peak, ifs_subrois["accs"].filter(regex="IFS$"), .33, (.32, .5, 7), False) f.text(.01, .97, "A", size=12) f.text(.01, .42, "B", size=12) f.text(.48, .42, "C", size=12) f.text(.825, .42, "D", size=12) sns.despine() save_figure(f, "figure_4") dksort_figure_4() **Decoding results for the dimensions rules in prefrontal cortex.** (A) Time-resolved decoding results. Solid lines indicate mean decoding accuracy within each time bin, and error bands denote bootstrapped standard error across subjects. The horizontal dashed line shows the empirical measure of chance performance, derived from the permutation analysis. The vertical dashed line is placed at the onset of the first stimulus. (B) The height of each bar shows the percentile in the shuffled null distribution corresponding to the observed accuracy for each subject and region. The horizontal dashed line demarcates the criterion of significance at a (corrected) alpha = 0.05. Bars are sorted by height within region. (C and D) Decoding accuracies fit to data averaged across the timepoints at 3 and 5 seconds following stimulus onset. Points and error bars represent mean and bootstrapped standard error across subjects, respectively. ### Dimension rule decoding during cue period¶ Now test whether we can decode the dimension rules during the cue period. This is both a) possibly interesting on its own and b) serves as a control against certain interpretations of our main results. In : pfc_cue = dksort_decode(["IFS"], "dimension_cue")  Take the mean accuracy across subjects and look at the maximum over time. In : def dksort_cue_test(): accs = pfc_cue["accs"].mean(axis=0) peak_tp = np.argmax(accs) max_accs = pfc_cue["accs"].values[:, peak_tp] low, high = moss.ci(moss.bootstrap(max_accs), 95) peak_time = pfc_cue["accs"].columns[peak_tp] args = max_accs.mean(), peak_time, low, high print "Max accuracy: %.3f at %ds; 95%% CI: %.3f, %.3f" % args dksort_cue_test()  Max accuracy: 0.363 at 5s; 95% CI: 0.338, 0.389  Now look at the maximum t statistic across time. The p value is corrected for multiple comparisons, here just across time. In : pfc_cue["ttest"]  Out: mu sd t p tp ROI IFS 0.36 0.05 2.12 0.12 5 ### Decision rule decoding in PFC¶ Next we turn to the Decision Rules (i.e. whether "yes" means "same" or "different"). We just repeat the main steps from above. First load up all the data and do the heavy lifting on the analysis side of things. In : pfc_decision = dksort_decode(pfc_rois, "decision")  Print the results of the group t test against chance. In : pfc_decision["ttest"]  Out: mu sd t p tp ROI FPC 0.49 0.06 -0.97 1.00 1 IFG 0.48 0.06 -1.45 1.00 1 IFS 0.51 0.06 0.40 0.94 5 aIns 0.50 0.05 0.24 0.96 -1 aMFG 0.50 0.06 -0.32 1.00 5 pMFG 0.50 0.05 0.11 0.98 5 pSFS 0.50 0.03 -0.02 0.99 3 Now also look at whether any individual models were significant. In : pfc_decision["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())  Out: IFS aMFG pMFG FPC IFG aIns pSFS correction omni 1 0 0 0 0 1 0 time 3 1 1 0 0 4 1 In : other_decision = dksort_peak_decode("decision", other_rois, other_masks, pfc_decision["peak"]["IFS"])  Now we can plot the decision rule figure In : def dksort_figure_5(): f = plt.figure(figsize=(4.48, 5.5)) ax_ts = f.add_axes([.12, .53, .86, .45]) dksort_timecourse_figure(ax_ts, pfc_decision, (.43, .67, 7)) ax_sig = f.add_axes([.12, .08, .34, .34]) dksort_significance_figure(ax_sig, pfc_decision["signif"]) ax_pfc_peak = f.add_axes([.59, .08, .22, .34]) dksort_point_figure(ax_pfc_peak, pfc_decision["peak"], .5, (.43, .67, 7), True) ax_sub_peak = f.add_axes([.86, .08, .12, .34]) dksort_point_figure(ax_sub_peak, other_decision["accs"].filter(regex="IFS$"), .5, (.43, .67, 7), False)

f.text(.01, .97, "A", size=12)
f.text(.01, .42, "B", size=12)
f.text(.48, .42, "C", size=12)
f.text(.825, .42, "D", size=12)

sns.despine()
save_figure(f, "figure_5")

dksort_figure_5() Plot conventions are identical to those in Figure 4. The y-axis is scaled to span similar binomial probabilities as in Figure 4.

## Control signals in posterior neocortex¶

We now return to the dimension rules and build decoding models to predict these rules from patterns in posterior cortex ROIS (specifically IPS and OTC). We'll carry forward the IFS ROI as we are interested in how the profile of information in this region is similar or different to the other regions. Lacking a good name for this collection of ROIs, we'll call them net as they form a sort of goal-directed attention network.

### Spatiotemporal decoding in IPS and visual cortex¶

In :
net_dimension = dksort_decode(net_rois, "dimension")

In :
net_dimension["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())

Out:
IFS IPS OTC
correction
omni 11 15 15
time 13 15 15
In :
net_dimension["ttest"]

Out:
mu sd t p tp
ROI
IFS 0.46 0.06 7.95 0 3
IPS 0.56 0.10 8.87 0 3
OTC 0.61 0.07 14.76 0 5

### Spatial scale of context information¶

We are interested in dissociating the information within this network. First we inquire about its spatial resolution. We do so by computing the first-order spatial autocorrelation on the maps of learned coefficients. Basically, what happens is that you extract the coefficients and put them in the native image space. Then, shift the map one voxel along each axis and compute the correlation with itself (ignoring the voxels that get pushed outside of the mask). Because there are three dimension rules, the logistic regression model actually has three binary models with associated weights. So we end up with correlation value for each axis/model, which we average to obtain a score for each subject and ROI.

In :
def dksort_calculate_ar1():
# Set up the dataframe
roi_index = moss.product_index([net_rois, subjects], ["roi", "subj"])
ar1s = pd.Series(index=roi_index, name="ar1", dtype=float)
shifters = [(0, 0, 1), (0, 1, 0), (1, 0, 0)]

for roi in net_rois:
ds = mvpa.extract_group("dimension", roi, mask_name, frames, peak, "rt")
coefs = mvpa.model_coefs(ds, model, flat=False)

for subj, coef_s in zip(subjects, coefs):
ar1_vals = []

# Average over the three underlying models
for coef_c in coef_s.transpose(3, 0, 1, 2):

# Average over shifts in the x, y, and z direction
for shifter in shifters:
x, y, z = (locs + shifter).T
shifted = coef_c[x, y, z]
notnan = ~np.isnan(shifted)
r, p = stats.pearsonr(orig[notnan], shifted[notnan])
ar1_vals.append(r)

ar1s.loc[(roi, subj)] = np.mean(ar1_vals)

# Print descriptive statistics about the autocorrelation
for roi in net_rois:
m, sd = ar1s.loc[roi].describe()[["mean", "std"]]
print "%s Spatial AR1: M = %.3g (SD = %.2g)" % (roi, m, sd)

# Perform pairwise tests for each ROI combination
for roi_a, roi_b in itertools.combinations(net_rois, 2):
t, p = stats.ttest_rel(ar1s.loc[roi_a], ar1s.loc[roi_b])
dof = len(ar1s.loc[roi_a]) - 1
print "Test for %s vs. %s: t(%d) = %.3f; p = %.3g" % (roi_a, roi_b, dof, t, p)

ar1_df = ar1s.reset_index()
return ar1_df

ar1_df = dksort_calculate_ar1()

IFS Spatial AR1: M = 0.153 (SD = 0.035)
IPS Spatial AR1: M = 0.151 (SD = 0.043)
OTC Spatial AR1: M = 0.249 (SD = 0.034)
Test for IFS vs. IPS: t(14) = 0.133; p = 0.896
Test for IFS vs. OTC: t(14) = -8.835; p = 4.23e-07
Test for IPS vs. OTC: t(14) = -6.498; p = 1.41e-05


Perform an omnibus test over the ROI factor (this actually comes before the pairwise tests above)

In :
%%R -i ar1_df
m.ar1 = lmer(ar1 ~ roi + (roi | subj), ar1_df, REML=FALSE)
m.ar1.nest = lmer(ar1 ~ 1 + (roi | subj), ar1_df, REML=FALSE)
print(anova(m.ar1))
lr_test(m.ar1, m.ar1.nest, "ROI effect")

Analysis of Variance Table
Df   Sum Sq  Mean Sq F value
roi  2 0.036324 0.018162  42.595
Likelihood ratio test for ROI effect:
Chisq(2) = 28.48; p = 6.53e-07


### Temproal profile of context decoding¶

Next we ask about temporal dissociations between the regions. This is hypothesis-driven: we expect that the IFS and IPS will have information about the attended dimensions earlier than the OTC. We're going to operationalize "time of information" with "time of peak information". There are other measures that might be insightful in different ways, but I think this is the best measurement given the type of data we have.

For a potentially more precise test, we are going to use cubic spline interpolation to upsample the original timeserieses to 500 ms bins. Then, we'll repeat the decoding analysis on these new timepoints. Finally, we'll fit gamma probability density functions to the high-res timecourse data and use those to infer on the time of peak information.

In :
def dksort_upsample():
if dv is not None:
dv["up_timepoints"] = up_timepoints
roi_index = moss.product_index([net_rois, subjects], ["ROI", "subj"])
columns = pd.Series(up_timepoints, name="timepoints")
up_accs = pd.DataFrame(index=roi_index, columns=columns, dtype=float)
up_hrf_info = pd.DataFrame(index=roi_index, columns=["peak", "r2"], dtype=float)
up_hrfs = {}
bounds = dict(shape=(2, 8), coef=(0, None), loc=(-2.5, 2.5), scale=(0, 3))

for roi in net_rois:
ds = mvpa.extract_group("dimension", roi + "_500ms", mask, frames,
upsample=4, confounds="rt", dv=dv4)
roi_accs = mvpa.decode_group(ds, model, dv=dv)
up_accs.loc[roi, :] = roi_accs
hrfs = [moss.GammaHRF(loc=-1.5, bounds=bounds) for s in subjects]
mapper = map if dv is None else dv.map_sync
hrfs_ = mapper(lambda h, a: h.fit(up_timepoints, a), hrfs, roi_accs)
up_hrfs[roi] = hrfs_
for subj, hrf in zip(subjects, hrfs_):
up_hrf_info.loc[(roi, subj)] = hrf.peak_time_, hrf.fit_r2_

return up_hrfs, up_hrf_info, up_accs

up_hrfs, up_hrf_info, up_accs = dksort_upsample()


Report main information about the fits.

In :
def dksort_peak_report():
for roi, info in up_hrf_info.groupby(level="ROI"):
peak_time = info["peak"].mean()
low, high = moss.ci(moss.bootstrap(info["peak"]), 95)
print roi + "\n---"
print " Mean peak: %.2fs;  95%% CI: (%.2f, %.2f)" % (peak_time, low, high)
r2 = info["r2"].median()
low, high = moss.ci(moss.bootstrap(info["r2"], func=np.median), 95)
print " Median R^2: %.2f;  95%% CI: (%.2f, %.2f)\n" % (r2, low, high)

dksort_peak_report()

IFS
---
Mean peak: 3.71s;  95% CI: (3.47, 3.95)
Median R^2: 0.77;  95% CI: (0.72, 0.89)

IPS
---
Mean peak: 3.47s;  95% CI: (3.22, 3.74)
Median R^2: 0.89;  95% CI: (0.78, 0.92)

OTC
---
Mean peak: 4.42s;  95% CI: (4.19, 4.64)
Median R^2: 0.95;  95% CI: (0.92, 0.96)



Define a short function to compute them mean and 95% CI for the pairwise differences and apply it to eaach pair.

In :
def up_peak_test(info, roi_a, roi_b):
peaks = np.array(info["peak"].unstack(level="ROI")[[roi_a, roi_b]])
diff = np.diff(peaks, axis=1)
mean_diff = np.mean(diff)
diff_ci = moss.ci(moss.bootstrap(diff), 95)
t, p = stats.ttest_rel(peaks[:, 0], peaks[:, 1])
print "%s - %s difference:" % (roi_a, roi_b)
print "  mean = %.2fs;" % mean_diff,
print "95%% CI: %.2f, %.2f" % tuple(diff_ci)
print "  t(%d) = %.3f; p = %.3g" % (len(diff) - 1, t, p)
print "  %d subjects with positive difference" % (diff > 0).sum()

In :
up_peak_test(up_hrf_info, "IFS", "OTC")

IFS - OTC difference:
mean = 0.71s; 95% CI: 0.41, 0.99
t(14) = -4.670; p = 0.000361
13 subjects with positive difference

In :
up_peak_test(up_hrf_info, "IPS", "OTC")

IPS - OTC difference:
mean = 0.95s; 95% CI: 0.77, 1.11
t(14) = -10.358; p = 6.03e-08
15 subjects with positive difference

In :
up_peak_test(up_hrf_info, "IFS", "IPS")

IFS - IPS difference:
mean = -0.24s; 95% CI: -0.53, 0.04
t(14) = 1.577; p = 0.137
3 subjects with positive difference


#### Now we can plot the figure corresponding to these analyses¶

In :
def dksort_figure_6():
f = plt.figure(figsize=(4.48, 5.5))

net_colors = [roi_colors[roi] for roi in net_rois]

# Set up the axes
ax_ts = f.add_axes([.12, .53, .855, .44])
ax_peak = f.add_axes([.12, .07, .16, .35])
ax_ar1 = f.add_axes([.37, .29, .61, .13])
ax_time = f.add_axes([.37, .07, .61, .13])

# Time-resolved decoding accuracy figure
# Plot the original-resolution accuracy as point estimate and CI
dksort_timecourse_figure(ax_ts, net_dimension, (.3, .66, 5),
"ci_bars", False, {"linewidth": 2}, ms=5)

# Plot the Gamma HRF modeled-accuracy as a smooth curve
xx = np.linspace(-2, 10, 500)
net_rois_r = list(reversed(net_rois))
fits = [[hrf.predict(xx) for hrf in up_hrfs[roi]] for roi in net_rois_r]
fits = np.transpose(fits, (1, 2, 0))
sns.tsplot(fits, time=xx, condition=net_rois_r, err_style=None, linewidth=1.75,
color=reversed(net_colors), ax=ax_ts)

# Remove the old chance line that doesn't span the whole plot
ax_ts.lines.remove()
ax_ts.plot([-2, 10], [1 / 3, 1 / 3], ls="--", c="k")

# Peak average decoding accuracy
dksort_point_figure(ax_peak, net_dimension["peak"], .33, (.3, .66, 5), True)

# Spatial autocorrelation of model parameters
box_colors = [sns.set_hls_values(c, l=.5, s=.3) for c in net_colors]
ar1_vals = [df.ar1.values for _, df in ar1_df.groupby("roi")]
sns.boxplot(ar1_vals, vert=False, color=box_colors, linewidth=1, widths=.6, ax=ax_ar1)
ax_ar1.set_yticklabels(net_rois)
ax_ar1.set_xlim(.05, .35)
ax_ar1.set_xticks([.1, .2, .3])
ax_ar1.xaxis.grid(True)
ax_ar1.set_xlabel("Spatial autocorrelation of feature weights (r)")

# Pairwise differences in peak decoding time
time_vals = [up_hrf_info.loc[roi, "peak"] - up_hrf_info.loc["OTC", "peak"] for roi in ["IFS", "IPS"]]
sns.boxplot(time_vals, vert=False, ax=ax_time, linewidth=1, widths=.5, color=box_colors[:2])
ax_time.set_xlim(-2.2, .7)
ax_time.xaxis.grid(True)
ax_time.set_yticklabels(net_rois[:2])
ax_time.axvline(0, ls=":", c="#444444")
ax_time.set_xlabel("Time of peak decoding relative to OTC (s)")

# Panel labels
text_size = 12
f.text(.01, .97, "A", size=text_size)
f.text(.01, .43, "B", size=text_size)
f.text(.3, .43, "C", size=text_size)
f.text(.3, .20, "D", size=text_size)

sns.despine()
save_figure(f, "figure_6")

dksort_figure_6()

/Users/mwaskom/anaconda/envs/dksort/lib/python2.7/site-packages/matplotlib/axes.py:4747: UserWarning: No labeled objects found. Use label='...' kwarg on individual plots.
warnings.warn("No labeled objects found. " **Decoding results for the dimension rules in posterior neocortex.** (A) Points and error bars show the mean and bootstrapped standard error for decoding accuracy in the original time bins. The solid traces show the gamma PDF models used to derive temporal information, averaged across subjects. Plot conventions are otherwise as in Figure 4A. (B) Plot conventions are as in Figure 4C. (C) Boxplots showing the distribution of model coefficient autocorrelation across subjects, sorted by region. (D) Boxplots showing the distribution of relative differences in the time of peak decoding accuracy for the IFS and IPS models relative to the OTC models. Negative numbers indicate later peaks in the OTC.

### Frontoparietal influence on visual selection¶

The next set of analyses focus on the question of influence, or at least interactions, within this network. Our question is, when the representation is stronger in frontoparietal areas (as measured by higher probabilistic prediction for the correct class), is that trial more likely to be classified correctly in OTC? I've actually done this with both binary and continuous DVs; the latter is actually a stronger effect, but the former provides a more interpreteable plot, with probabiltiy of correct prediction as the DV.

In :
def dksort_model_logits():
logits = dict()
preds = dict()
bins = np.arange(-14.5, 14.5, 1)
for roi in net_rois:
ds = mvpa.extract_group("dimension", roi, mask, frames, peak, "rt", dv=dv4)
logits_ = mvpa.decode_group(ds, model, logits=True, trialwise=True, dv=dv)
logits[roi] = np.concatenate(logits_)
preds_ = mvpa.decode_group(ds, model, trialwise=True, dv=dv)
preds[roi] = np.concatenate(preds_)

logit_df = behav_df[["subj"]]
for roi in net_rois:
logit_df[roi] = logits[roi]
logit_df[roi + "_acc"] = preds[roi]
logit_df[roi + "_bin"] = bins[np.digitize(logits[roi], bins)] - 1

return logit_df.reset_index(drop=True)
logit_df = dksort_model_logits()


Set up the main mixed model and the nested models for likelihood ratio tests.

In :
%%R -i logit_df
m.logits = lmer(OTC_acc ~ IFS + IPS + (IFS + IPS | subj), logit_df, family=binomial)
m.logits.noifs = lmer(OTC_acc ~ IPS + (IFS + IPS | subj), logit_df, family=binomial)
m.logits.noips = lmer(OTC_acc ~ IFS + (IFS + IPS | subj), logit_df, family=binomial)


Report on the model.

In :
%%R
print(m.logits, corr=FALSE)
lr_test(m.logits, m.logits.noifs, "IFS effect")
lr_test(m.logits, m.logits.noips, "IPS effect")

Generalized linear mixed model fit by the Laplace approximation
Formula: OTC_acc ~ IFS + IPS + (IFS + IPS | subj)
Data: logit_df
AIC  BIC logLik deviance
5062 5118  -2522     5044
Random effects:
Groups Name        Variance   Std.Dev.  Corr
subj   (Intercept) 4.7257e-02 0.2173867
IFS         3.0514e-05 0.0055239 -1.000
IPS         2.4432e-03 0.0494288 -0.353  0.353
Number of obs: 3961, groups: subj, 15

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.51027    0.06630   7.697 1.40e-14 ***
IFS          0.08242    0.01289   6.396 1.59e-10 ***
IPS          0.17603    0.01924   9.151  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Likelihood ratio test for IFS effect:
Chisq(1) = 24.49; p = 7.47e-07
Likelihood ratio test for IPS effect:
Chisq(1) = 29.10; p = 6.88e-08


It looks like the regression coefficient for the IPS is larger than the IFS. Is this statistically significant?

In :
%%R
C = matrix(c(0, -1, 1), 1, 3)
rownames(C) = "IPS-IFS"
print(summary(glht(m.logits, C)))

	 Simultaneous Tests for General Linear Hypotheses

Fit: glmer(formula = OTC_acc ~ IFS + IPS + (IFS + IPS | subj), data = logit_df,
family = binomial)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
IPS-IFS == 0  0.09361    0.02453   3.816 0.000136 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)



Now collect and add the motion estimates. We'll take the average relative displacement metric on the two TRs that the BOLD data were taking from on each trial.

In :
def dksort_collect_motion_info():
all_motion = []
motion_template = op.join(anal_dir, "dksort/%s/preproc/run_%d/realignment_params.csv")
for subj in subjects:
subj_motion = []
for run in range(1, 5):
motion = pd.read_csv(motion_template % (subj, run))["displace_rel"]
stim_onsets = behav_df.loc[(behav_df.subj == subj) & (behav_df.run == run), "stim_time"]
stim_onsets = stim_onsets.values.astype(int) / 2
first_tr = stim_onsets + 2
second_tr = stim_onsets + 3
run_motion = np.mean([motion[first_tr], motion[second_tr]], axis=0)
subj_motion.append(run_motion)
subj_motion = np.concatenate(subj_motion)
all_motion.append(subj_motion)
logit_df["motion"] = stats.zscore(np.concatenate(all_motion))
dksort_collect_motion_info()

In :
%%R -i logit_df
m.logits.motion = lmer(OTC_acc ~ IFS + IPS + motion + (IFS + IPS + motion | subj), logit_df, family=binomial)
print(m.logits.motion)

Generalized linear mixed model fit by the Laplace approximation
Formula: OTC_acc ~ IFS + IPS + motion + (IFS + IPS + motion | subj)
Data: logit_df
AIC  BIC logLik deviance
5061 5149  -2516     5033
Random effects:
Groups Name        Variance   Std.Dev.  Corr
subj   (Intercept) 4.7142e-02 0.2171227
IFS         1.4172e-05 0.0037646 -1.000
IPS         2.9784e-03 0.0545748 -0.308  0.308
motion      2.9775e-02 0.1725535 -0.464  0.464 -0.700
Number of obs: 3961, groups: subj, 15

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.50026    0.06663   7.508 6.02e-14 ***
IFS          0.08159    0.01291   6.320 2.61e-10 ***
IPS          0.17660    0.02019   8.747  < 2e-16 ***
motion      -0.07409    0.05872  -1.262    0.207
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) IFS    IPS
IFS     0.070
IPS    -0.167 -0.136
motion -0.273  0.056 -0.364


Plot a figure corresponding to these analyses

In :
def dksort_figure_7():

# Fit a logistic regression for each subject and save predictions
xx = np.linspace(-10, 6, 100)
pred_rois = ["IFS", "IPS"]
models = {roi: np.empty((subjects.size, xx.size)) for roi in pred_rois}
xlim = [-6, 4]

models = dict()
for roi in pred_rois:
fit = sm.GLM(logit_df.OTC_acc.values, x_fit, family=sm.families.Binomial()).fit()
models[roi] = fit.predict(x_pred)

f = plt.figure(figsize=(3.34, 5.5))
axr = f.add_axes([.15, .38, .8, .58])

# Plot the data for each ROI
for i, roi in enumerate(pred_rois):
color = roi_colors[roi]

histcolor = sns.set_hls_values(color, l=.5, s=.3)
sns.tsplot(models[roi], time=xx, color=color, label=roi,
linewidth=1.25, err_style=None, ax=axr)
roi_pivot = pd.pivot_table(logit_df, "OTC_acc", "subj", roi + "_bin").values
bins = np.sort(logit_df[roi + "_bin"].unique())
sns.tsplot(roi_pivot, time=bins, err_style="ci_bars", color=color, interpolate=False,
estimator=stats.nanmean, err_kws={"linewidth": 2}, markersize=6, ax=axr)
axh = f.add_axes([.15, .23 - i * .15, .8, .12])
bins = np.linspace(*xlim, num=xlim - xlim + 1)
axh.hist(logit_df[roi], bins, rwidth=.95,
color=histcolor, linewidth=0, label=roi)
axh.set_yticks(np.linspace(0, 900, 3))
axh.set_ylabel("Trials")
if i:
axh.set_xlabel("Classifier evidence for target class")
else:
axh.set_xticklabels([])
#axh.legend(loc="upper left")

# Set the regession axis properties
axr.set_xlim(*xlim)
axr.set_ylim(.25, .85)
axr.set_xticklabels([])
axr.axhline(.33, c="k", ls="--")
axr.set_ylabel("OTC classifier accuracy")

# Label the traces and histograms
axr.text(-3, .6, "IFS", color=roi_colors["IFS"], ha="center")
axr.text(-1.5, .48, "IPS", color=roi_colors["IPS"], ha="center")
f.axes.text(-4, 600, "IFS evidence", color=roi_colors["IFS"], ha="center")
f.axes.text(-4, 600, "IPS evidence", color=roi_colors["IPS"], ha="center")

# Label the panels
textsize = 11
f.text(.02, .96, "A", size=textsize)
f.text(.02, .35, "B", size=textsize)
f.text(.02, .20, "C", size=textsize)

sns.despine()
save_figure(f, "figure_7")

dksort_figure_7()