binned_spiking_activity_human_epilepsy/python_code/analysis_PlosCB/hierarchical_model.py

from pymer4.models import Lmer


def hierarchical_model_regions(dataframe):
    # overall p-value
    df_reduced = dataframe[['logepsilon', 'patient', 'binning_focus', 'region']]
    model = Lmer('logepsilon ~ binning_focus*region + (1|patient)',
                     data=df_reduced)

    model.fit()
    anova_result = model.anova(force_orthogonal=True)
    global_pvals = anova_result['P-val']
    # separate p-values
    pvals_regions = dict()
    for reg in dataframe.region.unique():
        print(reg)
        df_region = dataframe[dataframe['region'] == reg]
        df_reduced = df_region[['logepsilon', 'patient', 'binning_focus']]
        # for formulas, check
        # https://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet
        model = Lmer('logepsilon ~ binning_focus + (1|patient)',
                     data=df_reduced)
        model.fit()
        anova_result = model.anova(force_orthogonal=True)
        region_pval = anova_result['P-val'].item()
        pvals_regions[reg] = region_pval

    return global_pvals, pvals_regions


def hierarchical_model_hemispheres(dataframe):
    df_reduced = dataframe[['logepsilon', 'patient', 'binning_focus']]
    model = Lmer('logepsilon ~ binning_focus + (1|patient)',
                 data=df_reduced)
    model.fit()
    anova_result = model.anova(force_orthogonal=True)
    pval = anova_result['P-val'].item()
    return pval