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 [1]:
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 [2]:
import lyman
from lyman import mvpa, evoked
import seaborn as sns
import moss

Set up display options for figures and tables

In [3]:
%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 [4]:
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 [5]:
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 [6]:
%load_ext rmagic
%R library(lme4)
%R library(multcomp)
Loading required package: Matrix
Loading required package: lattice

Attaching package: ‘lme4’

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

    AIC, BIC

Loading required package: mvtnorm
Loading required package: survival
Loading required package: splines

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 [7]:
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 [8]:
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 [9]:
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 [10]:
subjects = pd.Series(lyman.determine_subjects(), name="subj")

Get root paths for lyman outputs

In [11]:
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 [12]:
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 [13]:
%%R
lr_test = function(m1, m2, name){
    out = anova(m1, m2)
    chi2 = out$Chisq[2]
    dof = out$"Chi Df"[2]
    p = out$"Pr(>Chisq)"[2]
    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 [14]:
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 [15]:
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 [16]:
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 [17]:
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 [18]:
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()[0])

    # 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[0])
            signif.loc[("omni", subj), roi] = stats.percentileofscore(null.max(axis=1), score[0])    

    return dict(rois=rois, accs=accs, chance=chance, peak=peak_df,
                null=null, ttest=ttest, signif=signif)
In [19]:
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 [20]:
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[1] - .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 [21]:
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[1]) + .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 [22]:
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 [23]:
%%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 [24]:
%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 [25]:
%%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 [26]:
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 [27]:
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 [28]:
%%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 [29]:
%%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 [30]:
%%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 [31]:
%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 [32]:
%%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 [33]:
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 [34]:
%%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 [35]:
%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 [36]:
%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 [37]:
%%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 [38]:
%%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)
lr_test(m.shift.add, m.shift.int, "interaction")
Likelihood ratio test for interaction:
  Chisq(1) = 0.61; p = 0.436
In [39]:
%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 [40]:
%%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)
lr_test(m.shift.add, m.shift.nodim, "dimension rule")
lr_test(m.shift.add, m.shift.nodec, "decision rule")
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 [41]:
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[0]
        y += text_offset[1]
        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 [42]:
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 [43]:
def dksort_roi_sizes():
    roi_sizes = pd.DataFrame(columns=all_rois, index=subjects, dtype=float)
    mask_template = op.join(data_dir, "%s/masks/yeo17_%s.nii.gz")
    for roi in roi_sizes:
        for subj in subjects:
            vox_count = nib.load(mask_template % (subj, roi.lower())).get_data().sum()
            roi_sizes.loc[subj, roi] = vox_count
    return roi_sizes.describe().loc[["mean", "std"]].astype(int)
dksort_roi_sizes()
Out[43]:
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 [44]:
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[44]:
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 [45]:
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 [46]:
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 [47]:
pfc_dimension["ttest"]
Out[47]:
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 [48]:
pfc_dimension["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())
Out[48]:
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 [49]:
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 [50]:
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)
    for roi, mask in zip(rois, masks):
        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 [51]:
ifs_subrois = dksort_peak_decode("dimension", other_rois, other_masks, pfc_dimension["peak"]["IFS"])
In [52]:
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 [53]:
ifs_subrois["accs"].describe().T
Out[53]:
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 [54]:
ifs_subrois["accs"].apply(lambda d: pd.Series(moss.ci(moss.bootstrap(d), 95), index=[2.5, 97.5]), axis=0)
Out[54]:
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 [55]:
ifs_subrois["ttest"].T
Out[55]:
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 [56]:
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 [57]:
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 [58]:
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 [59]:
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 [60]:
pfc_cue = dksort_decode(["IFS"], "dimension_cue")

Take the mean accuracy across subjects and look at the maximum over time.

In [61]:
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 [62]:
pfc_cue["ttest"]
Out[62]:
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 [63]:
pfc_decision = dksort_decode(pfc_rois, "decision")

Print the results of the group t test against chance.

In [64]:
pfc_decision["ttest"]
Out[64]:
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 [65]:
pfc_decision["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())
Out[65]:
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 [66]:
other_decision = dksort_peak_decode("decision", other_rois, other_masks, pfc_decision["peak"]["IFS"])

Now we can plot the decision rule figure

In [67]:
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 [68]:
net_dimension = dksort_decode(net_rois, "dimension")
In [69]:
net_dimension["signif"].groupby(level="correction").apply(lambda x: (x > 95).sum())
Out[69]:
IFS IPS OTC
correction
omni 11 15 15
time 13 15 15
In [70]:
net_dimension["ttest"]
Out[70]:
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 [71]:
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:
        mask_name = "yeo17_" + roi.lower()
        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):
                mask = ~np.isnan(coef_c)
                orig = coef_c[mask]
                locs = np.argwhere(mask)
                
                # 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 [72]:
%%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 [73]:
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:
        mask = "yeo17_" + roi.lower()
        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 [74]:
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 [75]:
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 [76]:
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 [77]:
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 [78]:
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 [79]:
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[0].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 [80]:
def dksort_model_logits():
    logits = dict()
    preds = dict()
    bins = np.arange(-14.5, 14.5, 1)
    for roi in net_rois:
        mask = "yeo17_" + roi.lower()
        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 [81]:
%%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 [82]:
%%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 [83]:
%%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 [84]:
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 [85]:
%%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 [86]:
def dksort_figure_7():
    
    # Fit a logistic regression for each subject and save predictions
    xx = np.linspace(-10, 6, 100)
    x_pred = sm.add_constant(xx, prepend=True)
    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:
        x_fit = sm.add_constant(logit_df[roi].values, prepend=True)
        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[1] - xlim[0] + 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[1].text(-4, 600, "IFS evidence", color=roi_colors["IFS"], ha="center")
    f.axes[2].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()