Syllable statistical testing

For comparing significant syllables between experimental conditions, we typically run a z-test on bootstrapped syllable usage distributions (one test per syllable), then correct for multiple comparisons using the Benjamini-Hochberg procedure at a false discovery rate of 5-10%. This procedure is best carried out in a Jupyter notebook.

Step 1. Load data

from moseq2_viz.model.util import parse_model_results
from moseq2_viz.util import parse_index
import numpy as np

#point paths to model and index files
model_file = '/path/to/model/file/file.p'
index_file = '/path/to/index/file/moseq2-index.yaml' #generate using "moseq2-viz generate-index"

#parse index and model results
index, sorted_index = parse_index(index_file)
model_results = parse_model_results(model_file)

#load data to pandas data frame
df = scalars_to_dataframe(sorted_index, include_model = model_file)

Step 2. Organize data

#indicate which subject names belong to which experimental group
group1 = ['C57-1_control_2019-01-01',
'C57-2_control_2019-01-01',
'C57-3_control_2019-01-01',
'C57-4_control_2019-01-01',
'C57-5_control_2019-01-01']

group2 = ['C57-1_experimental_2019-01-01',
'C57-2_experimental_2019-01-01',
'C57-3_experimental_2019-01-01',
'C57-4_experimental_2019-01-01',
'C57-5_experimental_2019-01-01']

#update df to include group labels for associated subject names
for _id in group1:
    df.loc[df.SubjectName==_id, 'group'] = 'group1'
for _id in group2:
    df.loc[df.SubjectName==_id, 'group'] = 'group2'

Step 3. Compute label usages

def gi(x):
    return ~np.isnan(x).any(tuple(range(1,x.ndim)))
def good(x):
    return x[gi(x)]
def calc_label_usage(labels, n_states=None):
    usage = np.bincount(np.nan_to_num(good(labels)).astype('int32'), minlength=n_states)
    usage = usage / float(usage.sum())
    return usage

#N=number of syllables to analyze — typically the number of syllables that account for 90% of all behavior
group1_usages = []
for i in df.SubjectName[df.group == 'group1'].unique():
    working_labels = df[df.SubjectName == i]['model_label (sort=frames)']
    group1_usages.append(calc_label_usage(working_labels[working_labels>-5].values)[:N])
    
group2_usages = []
for i in df.SubjectName[df.group == 'group2'].unique():
    working_labels = df[df.SubjectName == i]['model_label (sort=frames)']
    group2.append(calc_label_usage(working_labels[working_labels>-5].values)[:N])

Step 4. Bootstrap usages

def bootstrap_group(lst, rng):
    return list(rng.choice(len(lst),len(lst),replace=True))

def bootstrap_me(usages, iters=1000):
    bootstrap_mean_usages = []
    for i in range(iters):
        rng = np.random.RandomState(seed=i)
        temp = []
        boot_mice = bootstrap_group(usages, rng)
        for mouse in boot_mice:
            temp.append(usages[mouse])
        temp = np.asarray(temp)
        bootstrap_mean_usages.append(np.nanmean(temp,axis=0))
        
    return bootstrap_mean_usages

group1_boots = bootstrap_me(group1_usages)
group2_boots = bootstrap_me(group2_usages)

Step 5. Z-test

from scipy import stats

def ztest(d1, d2, mu1=None, mu2=None):
    mu1 = d1.mean() if mu1 is None else mu1
    mu2 = d2.mean() if mu2 is None else mu2
    std1, std2 = d1.std(), d2.std()

    std = np.sqrt(std1**2 + std2**2)
    return np.minimum(1.,2*stats.norm.cdf(-np.abs(mu1 - mu2)/std))

# do a ztest on the bootstrap distributions of your 2 conditions
pvals_ztest_boots = np.array([ztest(group1_boots[:,i], group2_boots[:,i]) for i in range(group1_boots.shape[1])])

Step 6. Multiple comparisons p-value correction (need statsmodels library)

from statsmodels.stats.multitest import multipletests

# significant syllables (relabeled by time used)
np.where(multipletests(pvals_ztest_boots, alpha=0.10, method='fdr_bh')[0])[0]