Python MultiComparison.allpairtest方法代码示例

# 需要导入模块: from statsmodels.stats.multicomp import MultiComparison [as 别名]
# 或者: from statsmodels.stats.multicomp.MultiComparison import allpairtest [as 别名]
def position_stats(df, name_mapping=None):

    # print '### position stats'
    from statsmodels.stats.weightstats import ztest
    from functools32 import partial, wraps
    POS = df.position.unique()
    POS.sort()
    model = 'value ~ group'
    allpvals = None
    header = None
    DF = None

    ttest_log_wrap = wraps(
        partial(ttest_ind_log, equal_var=False))(ttest_ind_log)
    ttest_ind_nev = wraps(
        partial(stats.ttest_ind, equal_var=False))(stats.ttest_ind)
    mwu_test = wraps(partial(stats.mannwhitneyu, use_continuity=False))(
        stats.mannwhitneyu)

    bootstrap_sample_num = 1000
    # print df

    stats_test = ttest_ind_nev
    GROUPS = df.group.unique()
    # GROUPS = [0,3]

    for pos in POS:
        # print pos
        data = df[df.position == pos]
        data = data.groupby(['sid']).mean()
        data = resample_data(data, num_sample_per_pos=BOOTSTRAP_NUM)
        # print data
        # print data.group.unique()
        # data = df[(df.group == 0) | (df.group == 3)]
        # print data
        # sys.exit()

        #cross = smf.ols(model, data=data).fit()
        #anova = sm.stats.anova_lm(cross, type=1)
        # print data.group

        mcp = MultiComparison(data.value, data.group.astype(int))

        rtp = mcp.allpairtest(stats_test, method='bonf')
        mheader = []
        for itest in rtp[2]:
            name1 = itest[0]
            name2 = itest[1]
            if name_mapping is not None:
                name1 = name_mapping[str(name1)]
                name2 = name_mapping[str(name2)]

            mheader.append("{} - {}".format(name1, name2))

        if not header or len(mheader) > len(header):
            header = mheader

        # get the uncorrecte pvals
        pvals = rtp[1][0][:, 1]

        ndf = pd.DataFrame(data=[pvals], columns=mheader)
        if allpvals is None:
            allpvals = ndf
        else:
            allpvals = pd.concat([allpvals, ndf])

    # return allpvals
    # corr_pvals = allpvals
    # print allpvals
    # return allpvals

    flatten = allpvals.values.ravel()
    flatten = flatten * 2
    mcpres = multipletests(flatten, alpha=0.05, method='bonf')
    # print mcpres
    corr_pvals = np.array(mcpres[1])
    # print corr_pvals
    corr_pvals = np.reshape(corr_pvals, (len(POS), -1))

    # print corr_pvals,corr_pvals.shape,header
    data = pd.DataFrame(data=corr_pvals, columns=header)
    data = data[data.columns[:3]]
    return data

# 需要导入模块: from statsmodels.stats.multicomp import MultiComparison [as 别名]
# 或者: from statsmodels.stats.multicomp.MultiComparison import allpairtest [as 别名]
plt.xlim(*xlim)
pair_labels = mod.groupsunique[np.column_stack(res2[1][0])]
plt.xticks([0,1,2], pair_labels)
plt.title('Multiple Comparison of Means - Tukey HSD, FWER=0.05' +
          '\n Pairwise Mean Differences')          

# Save to outfile
outFile = 'MultComp.png'
plt.savefig('MultComp.png', dpi=200)
print 'Figure written to {0}'.format(outFile)

plt.show()

# Instead of the Tukey's test, we can do pairwise t-test
# First, with the "Holm" correction
rtp = mod.allpairtest(stats.ttest_rel, method='Holm')
print rtp[0]

# and then with the Bonferroni correction
print mod.allpairtest(stats.ttest_rel, method='b')[0]

# Done this way, the variance is calculated at each comparison.
# If you want the joint variance across all samples, you have to 
# use a few tricks:(http://jpktd.blogspot.co.at/2013/03/multiple-comparison-and-tukey-hsd-or_25.html)
res2 = pairwise_tukeyhsd(dta2['StressReduction'], dta2['Treatment'])
studentized_mean = res2[1][2]
studentized_variance = res2[1][3]

t_stat = (studentized_mean / studentized_variance) / np.sqrt(2)
dof = len(dta2) - len(mod.groupsunique)
my_pvalues = stats.t.sf(np.abs(t_stat), dof) * 2  # two-sided

# 需要导入模块: from statsmodels.stats.multicomp import MultiComparison [as 别名]
# 或者: from statsmodels.stats.multicomp.MultiComparison import allpairtest [as 别名]
def main():
    # Note: the statsmodels module is required here.
    from statsmodels.stats.multicomp import (pairwise_tukeyhsd,
                                             MultiComparison)
    from statsmodels.formula.api import ols
    from statsmodels.stats.anova import anova_lm
    
    # Set up the data, as a structured array.
    # The first and last field are 32-bit intergers; the second field is an
    # 8-byte string. Note that here we can also give names to the individual
    # fields!
    dta2 = np.rec.array([
    (  1,   'mental',  2 ),
    (  2,   'mental',  2 ),
    (  3,   'mental',  3 ),
    (  4,   'mental',  4 ),
    (  5,   'mental',  4 ),
    (  6,   'mental',  5 ),
    (  7,   'mental',  3 ),
    (  8,   'mental',  4 ),
    (  9,   'mental',  4 ),
    ( 10,   'mental',  4 ),
    ( 11, 'physical',  4 ),
    ( 12, 'physical',  4 ),
    ( 13, 'physical',  3 ),
    ( 14, 'physical',  5 ),
    ( 15, 'physical',  4 ),
    ( 16, 'physical',  1 ),
    ( 17, 'physical',  1 ),
    ( 18, 'physical',  2 ),
    ( 19, 'physical',  3 ),
    ( 20, 'physical',  3 ),
    ( 21,  'medical',  1 ),
    ( 22,  'medical',  2 ),
    ( 23,  'medical',  2 ),
    ( 24,  'medical',  2 ),
    ( 25,  'medical',  3 ),
    ( 26,  'medical',  2 ),
    ( 27,  'medical',  3 ),
    ( 28,  'medical',  1 ),
    ( 29,  'medical',  3 ),
    ( 30,  'medical',  1 )], dtype=[('idx', '<i4'),
                                    ('Treatment', '|S8'),
                                    ('StressReduction', '<i4')])
    
    # First, do an one-way ANOVA
    df = pd.DataFrame(dta2)
    model = ols('StressReduction ~ C(Treatment)',df).fit()
    
    anovaResults =  anova_lm(model)
    print(anovaResults)
    if anovaResults['PR(>F)'][0] < 0.05:
        print('One of the groups is different.')
    
    #Then, do the multiple testing
    mod = MultiComparison(dta2['StressReduction'], dta2['Treatment'])
    print((mod.tukeyhsd().summary()))
    
    # The following code produces the same printout
    res2 = pairwise_tukeyhsd(dta2['StressReduction'], dta2['Treatment'])
    #print res2[0]
    
    # Show the group names
    print((mod.groupsunique))
    
    # Generate a print
    import matplotlib.pyplot as plt
    xvals = np.arange(3)
    plt.plot(xvals, res2.meandiffs, 'o')
    #plt.errorbar(xvals, res2.meandiffs, yerr=np.abs(res2[1][4].T-res2[1][2]), ls='o')
    errors = np.ravel(np.diff(res2.confint)/2)
    plt.errorbar(xvals, res2.meandiffs, yerr=errors, ls='o')
    xlim = -0.5, 2.5
    plt.hlines(0, *xlim)
    plt.xlim(*xlim)
    pair_labels = mod.groupsunique[np.column_stack(res2._multicomp.pairindices)]
    plt.xticks(xvals, pair_labels)
    plt.title('Multiple Comparison of Means - Tukey HSD, FWER=0.05' +
              '\n Pairwise Mean Differences')          
    
    # Save to outfile
    outFile = 'MultComp.png'
    plt.savefig('MultComp.png', dpi=200)
    print(('Figure written to {0}'.format(outFile)))
    
    plt.show()
    
    # Instead of the Tukey's test, we can do pairwise t-test
    # First, with the "Holm" correction
    rtp = mod.allpairtest(stats.ttest_rel, method='Holm')
    print((rtp[0]))
    
    # and then with the Bonferroni correction
    print((mod.allpairtest(stats.ttest_rel, method='b')[0]))
    
    # Done this way, the variance is calculated at each comparison.
    # If you want the joint variance across all samples, you have to 
    # use a few tricks:(http://jpktd.blogspot.co.at/2013/03/multiple-comparison-and-tukey-hsd-or_25.html)
    res2 = pairwise_tukeyhsd(dta2['StressReduction'], dta2['Treatment'])
    studentized_mean = res2.meandiffs
#.........这里部分代码省略.........