数据探索之假设检验

作者: 鱼心DrFish | 来源:发表于2017-06-16 17:13 被阅读633次

    推断统计学一般有两种方法,一是使用置信区间估算总体的参数,二是对总体参数的假设值进行决策。后者被称为假设检验,是我们这篇文章所要探讨的主题。

    本文首先会介绍假设检验的逻辑,然后利用BRFSS数据探讨富人是否更瘦这个问题,除了使用t检验方法外,还采用了蒙特卡洛模拟的方法来实现假设检验。

    假设检验的逻辑

    假设检验的逻辑有些类似于反证法,这里有原假设备择假设,它们之间是相互对立的,我们先在原假设成立的前提下进行计算,如果最终结果证明原假设不成立,我们就会接受备择假设,所以一般我们也会将想要证明的结论作为备择假设。

    比如,在讨论富人是否比普通人更瘦的问题上,原假设是:富人和普通人一般胖瘦;而备择假设是:富人比普通人更瘦。然后我们通过实际的样本数据,证明在原假设成立的前提下:如果得到现有样本结果(甚至更极端的情况)的可能性微乎其微,以至于我们需要拒绝原假设而接受备择假设;如果可能性在我们能够接受的范围内,我们则无法拒绝原假设。

    假设检验一般有下面几个步骤:

    • 设置原假设和备择假设。
    • 定义检验统计量。
    • 在原假设的基础上,通过计算或模拟得到检验统计量和P值。
    • 做判断:如果P值小于你能容忍的显著性水平时,则拒绝原假设。

    这里的检验统计量是由样本数据计算而来的,如果检验统计量的值出现在拒绝域中,则拒绝原假设。而一般我们会计算检验统计量出现在极端情况下的概率,也就是P值,并与事先设置的显著性水平比较,如果P值更小,则意味着检验统计量处在拒绝域,需要拒绝原假设。

    单个样本均值的假设检验

    在开始假设检验前,我们需要先导入相应的Python模块和BRFSS数据,其中采用了scipy的stats统计模块,它提供了假设检验计算P值的函数。

    import pandas as pd
    import numpy as np
    import matplotlib.pyplot as plt
    import seaborn as sns
    import brfss
    from scipy import stats  # 导入scipy的统计模块
    
    %config InlineBackend.figure_format = 'retina'
    
    df = brfss.ReadBrfss()  # 读取数据
    

    然后根据收入水平选取富人和普通人的bmi数据值。

    df2 = df[['bmi', 'income']].dropna()
    bmi = df2.bmi
    bmi_rich = df2[df2.income == 8].bmi  # 富人bmi数据
    bmi_ord = df2[df2.income != 8].bmi  # 普通人bmi数据
    

    我们需要判断的问题是:富人会更瘦吗?为了演示单个样本的情况,这里只使用富人bmi数据和所有人总体bmi均值做比较。并假设总体的bmi均值等于样本bmi均值,由于样本量有40万,根据中心极限定理可以做这样的近似。

    mu = np.mean(bmi) # 计算bmi均值
    mu

        28.18812531332513
    
    
    
    在只使用富人bmi数据的单样本情况下,原假设为:富人bmi均值等于总体bmi均值28.188;备择假设为:富人bmi均值小于总体bmi均值28.188,这是一个单边检验。这里采用t检验,所以检验统计量为t统计量。
    
    t检验一般用于变量满足或近似正态分布且总体标准差未知的情况下,特别是在小样本量的时候。但在大样本情况下,t检验和z检验是近似的,因此也可使用t检验。在scipy.stats模块中t检验的函数是`ttest_1samp()`,用于单个样本和总体均值比较的情况。其返回结果是t统计量和P值。
    
    
    
    >```python
    t_stats, p_value = stats.ttest_1samp(bmi_rich, mu)
    print("p value is %.10f" % (p_value/2))  # 除以2是因为我们使用的是单边检验
    
    p value is 0.0000000000
    

    以上使用了所有富人的bmi数据,得到P值在目前的精度下为0,所以我们拒绝原假设,接受备择假设,即得出结论:富人更瘦。如果我们只使用部分bmi数据,则会得到完全不一样的结果。

    t_stats, p_value = stats.ttest_1samp(bmi_rich[:500], mu) # 选择前500个数据进行检验
    print("p value is %.10f" % (p_value/2))

        p value is 0.2906513576
    
    
    
    >```python
    t_stats, p_value = stats.ttest_1samp(bmi_rich[:4000], mu)  # 选择前4000个数据进行检验
    print("p value is %.10f" % (p_value/2))
    
    p value is 0.0000660545
    

    当选择前500个数据进行t检验时,P值为0.29,而如果我们选择的显著性水平为0.1时,则无法拒绝原假设,但也不能证明原假设就成立。当选择前4000个数据进行t检验,P值为0.000066, 显著性水平为0.1时,则可以拒绝原假设。所以样本量越大,你得到的结果越明确。

    蒙特卡洛模拟之 bootstrap 方法

    以上直接采用了t检验,有时t检验的前提无法满足,我们还可以采用计算机模拟的方法,常用的是bootstrap方法。

    bootstrap是一种自助的方法,即在原有的样本中进行可重复的重新抽样,得到更多的样本。因为我们先假设富人和普通人的bmi值没有差别,所以从富人的样本数据中进行重新抽样,其实是模拟了普通人的样本分布。然后将原富人样本bmi均值与多个新样本的bmi均值比较,从而计算出极端情况下的概率,即P值。

    以下是实现bootstrap方法的代码。

    def bootstrap_replicate_1d(data, func):   # 进行一次重新抽样,并返回检验统计量
        return func(np.random.choice(data, size=len(data)))  
    
    def draw_bs_reps(data, func, size=1):
        bs_replicates = np.empty(size)  # 初始一个空数组
        for i in range(size):   # 进行多次重新抽样
            bs_replicates[i] = bootstrap_replicate_1d(data, func)  
        return bs_replicates  # 返回多次抽样的检验统计量数组
    
    
    def bootstrap_pvalue_1samp(data, pop_stats, func, size=1):
        sample_stats = func(data)  # 计算原有样本的检验统计量    
        translated_data = data - sample_stats + pop_stats  # 数据平移
        bs_replicates = draw_bs_reps(translated_data, func, size) # 重新抽样
        p = np.sum( bs_replicates < sample_stats) / size # 计算抽样统计量小于原有统计量的概率
        return p
    

    同样,我们使用三组不同大小的数据得到了不同的结果,在大样本的情况下P值很小,需要拒绝原假设。

    bootstrap_pvalue_1samp(bmi_rich, mu, np.mean, size=10000)

        0.0
    
    
    
    
    >```python
    bootstrap_pvalue_1samp(bmi_rich[:500], mu, np.mean, size=10000)
    
    0.71260000000000001
    

    bootstrap_pvalue_1samp(bmi_rich[:4000], mu, np.mean, size=10000)

        0.0001
    
    
    
    ### 两个独立样本均值的假设检验
    
    现在我们讨论两个独立样本的情况,即使用富人bmi数据和普通人bmi数据。原假设为:富人和普通人的bmi均值相同;备择假设为:富人bmi均值比普通人小,仍是单边检验。同样,我们采用t检验,使用scipy.stats中的`ttest_ind()`函数,计算两个独立样本t检验的P值。
    
    选择显著性水平为0.1时,当样本量只有500时,P值大于0.1,无法拒绝原假设。而当样本量为1000或全部时,P值小于0.1,拒绝原假设,证明了富人比普通人更瘦。
    
    
    >```python
    t_stats, p_value = stats.ttest_ind(bmi_rich, bmi_ord)
    print("p value is %.10f" % (p_value/2))
    
    p value is 0.0000000000
    

    t_stats, p_value = stats.ttest_ind(bmi_rich[:500], bmi_ord[:500])
    print("p value is %.10f" % (p_value/2))

        p value is 0.1219871502
    
    
    
    >```python
    t_stats, p_value = stats.ttest_ind(bmi_rich[:1000], bmi_ord[:1000])
    print("p value is %.10f" % (p_value/2))
    
    p value is 0.0000846839
    

    蒙特卡洛模拟之Permutation方法

    同样,在两个独立样本的情况下,我们也可以采用计算机模拟的方法,这里介绍另一种Permutation的方法。

    Permutation方法的逻辑是,在原假设成立的前提下,富人bmi和普通人bmi数据没有区别,所以可以将它们混合重新洗牌,然后重新分成两组数据,并计算这两组新数据均值的差值,重复多次这一操作。然后再将原有富人和普通人bmi均值之差与这些新组合的差值进行比较,计算极端情况下的概率,如果概率很小,则推翻了无差别的原假设。

    以下是实现Permutation方法的代码。

    def diff_of_means(data_1, data_2):  
        diff = np.mean(data_2) - np.mean(data_1)  # 计算两组数据均值的差
        return diff
    
    
    def permutation_sample(data1, data2):  # 产生新的分组数据
        data = np.concatenate((data1, data2))  # 合并两组数据
        permuted_data = np.random.permutation(data)  # 对合并后的数据进行重新排列
        perm_sample_1 = permuted_data[:len(data1)]  # 分成新的两组数据
        perm_sample_2 = permuted_data[len(data1):]
        return perm_sample_1, perm_sample_2
    
    
    def draw_perm_reps(data_1, data_2, func, size=1):  # 进行多次重新分组的操作
        perm_replicates = np.empty(size)
        for i in range(size):
            perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)
            perm_replicates[i] = func(perm_sample_1, perm_sample_2)
        return perm_replicates
    
    
    def permutation_pvalue(data_1, data_2, func, size=1):  # 计算P值
        empirical_test_stats = func(data_1, data_2)
        perm_replicates = draw_perm_reps(data_1, data_2, func, size)
        p = np.sum(perm_replicates > empirical_test_stats) / len(perm_replicates)
        return p
    
    
    

    同样,我们使用了大小不同的三组样本数据进行测试,得到与上述t分布一致的判断结果。

    permutation_pvalue(bmi_rich, bmi_ord, diff_of_means, size=10000)

        0.0
    
    
    
    
    >```python
    permutation_pvalue(bmi_rich[:500], bmi_ord[:500], diff_of_means, size=10000)
    
    0.1142
    

    permutation_pvalue(bmi_rich[:1000], bmi_ord[:1000], diff_of_means, size=10000)

        0.0001
    
    
    
    ### 小结
    
    至此,本系列的数据探索文章,使用美国行为风险因素监控BRFSS调研数据,围绕“富人会更瘦还是更胖” 这个问题,采用了多种统计学方法来寻找答案,它们分别是:
    
    * 计算bmi均值的差
    * 绘制统计图表
    * 计算Cohen's d
    * 计算bmi的置信区间
    * 假设检验方法
    
    由简到难,在解决问题的道路上,我们也将基础的统计方法在Python实现的情况下进行了整理和回顾,希望对大家有帮助。
    
    ---
    [本文代码github](https://github.com/fishstar/Exploratory-Data-Analysis)
    
    数据探索系列目录:
    1. [开篇:数据分析流程](http://www.jianshu.com/p/dffdaf11bd4c)
    2. [描述性统计分析](http://www.jianshu.com/p/8982ad63eb85)
    3. [统计分布](http://www.jianshu.com/p/8a0479f55b21)
    4. [参数估计](http://www.jianshu.com/p/7e556f17021a)
    5. 假设检验(本文)
    
    参考资料:
    * [《Think Stats 2》](http://greenteapress.com/wp/think-stats-2e/)
    * [《统计学》](https://book.douban.com/subject/4076214/),William Mendenhall著
    
    致谢:
    最后还是非常感谢解密大数据社群的小伙伴们给我的支持和鼓励,让我们一起成长。
    ![解密大数据社群公号](http:https://img.haomeiwen.com/i4420255/691af59970307eff.jpg?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

    相关文章

      网友评论

      • 2f29808a065a:“ 原假设是:富人和普通人一般胖瘦;而备择假设是:富人比普通人更瘦。”
        原假设 和 备择假设 应该是完备且对立。文中所陈述明显有问题。
        鱼心DrFish:@lzx_465c 是不太严谨

      本文标题:数据探索之假设检验

      本文链接:https://www.haomeiwen.com/subject/dwsmqxtx.html