方差分析介绍(结合COVID-19案例)

时间  2020-06-17
标签 方差 分析 介绍 结合 covid 案例 繁體版
原文   原文链接

做者|GUEST BLOG
编译|VK
来源|Analytics Vidhyapython

介绍

“事实是每一个人都相信的简单陈述。也就是事实是没有错的，除非它被人发现了错误。假设有一个没人愿意相信的建议，那么它要直到被发现有效的时候才能成为事实。” –爱德华·泰勒

咱们正在应对一场空前规模的流行病。全世界的研究人员都在疯狂地试图开发一种疫苗或COVID-19的治疗方法，而医生们正试图阻止这种流行病席卷整个世界。git

我最近有了一个想法，把个人统计知识应用到这些大量COVID数据中。github

考虑这样一个场景：医生有四种医疗方法来治疗病人。一旦咱们有了测试结果，用最少时间治愈病人的治疗会是最好的方法。api

但若是这些病人中的一些已经部分治愈，或者其余药物已经在治疗他们呢？app

为了做出一个有信心和可靠的决定，咱们须要证据来支持咱们的作法。这就是方差分析的概念发挥做用的地方。dom

在本文中，我将向你介绍方差分析测试及其用于作出更好决策的不一样类型。我将在Python中演示每种类型的ANOVA（方差分析）测试，以可视化它们并处理COVID-19数据。机器学习

注意：你必须了解统计学的基本知识才能理解这个主题。最好了解t检验和假设检验。函数

什么是方差分析测试(ANOVA)

方差分析，或称方差分析，能够看做是两组以上的t检验的推广。独立t检验用于比较两组之间的条件平均值。当咱们想比较两组以上患者的病情平均值时，使用方差分析。学习

方差分析测试模型中某个地方的平均值是否存在差别（测试是否存在总体效应），但它不能告诉咱们差别在哪里（若是存在）。为了找出两组之间的区别，咱们必须进行过后检验。测试

要执行任何测试，咱们首先须要定义原假设和替代假设：

  • 零假设–各组之间无显着差别
  • 替代假设–各组之间存在显着差别

基本上，方差分析是经过比较两种类型的变化来完成的，即样本均值之间的变化，以及每一个样本内部的变化。如下公式表示单向Anova测试统计数据。

ANOVA公式的结果，即F统计量（也称为F比率），容许对多组数据进行分析，以肯定样本之间和样本内部的可变性。

单向ANOVA的公式能够这样写：

当咱们绘制ANOVA表时，上面的全部组成部分均可以以下所示：

通常来讲，若是与F相关联的p值小于0.05，则将拒绝原假设并支持替代假设。若是原假设被拒绝，咱们能够得出结论，全部组的均值不相等。

注：若是被测组之间不存在真正的差别，也就是所谓的零假设，那么方差分析的F比统计结果将接近1。

ANOVA检验的假设

在进行方差分析以前，咱们须要作一些假设：

  1. 从因子水平定义的整体中独立且随机地得到观察结果
  2. 每一个因子水平的数据均呈正态分布
  3. 案例独立性：样本案例应相互独立
  4. 方差的同质性：同质性是指各组之间的方差应近似相等

方差同质性的假设能够用Levene检验或Brown-Forsythe检验来检验。分数分布的正态性能够用直方图、偏度和峰度值来检验，也能够用Shapiro-Wilk或Kolmogorov-Smirnov或Q-Q图来检验。独立性的假设能够根据研究设计来肯定。

值得注意的是，方差分析对于假设独立性的违规行为并不强大。这就是说，即便你违反了同质性或正态性的假设，你也能够进行测试并基本相信结果。

可是，若是违反了独立性假设，方差分析的结果是无效的。通常来讲，在违反同质性的状况下，若是具备相同大小的组，则分析被认为是可靠的。对于违反正态性的状况，若是样本量较大，继续进行方差分析一般是能够的。

方差分析检验类型

  1. 单向方差分析：单向方差分析只有一个自变量

    • 例如，能够按国家/地区评估日冕案例的差别，而且一个国家能够将2个，20个或更多不一样的类别进行比较

  2. 双向方差分析：双向方差分析（也称为因子方差分析）是指使用两个独立变量的方差分析

    • 扩展上面的示例，双向方差分析能够按年龄组（独立变量1）和性别（独立变量2）检查日冕病例（因变量）的差别。双向方差分析可用于检查两个独立变量之间的相互做用。相互做用代表，自变量的全部类别之间的差别不是统一的
    • 例如，老年组整体上可能比青年组具备更高的日冕病例，可是与欧洲国家相比，亚洲国家的差别可能更大（或更小）

  3. N向方差分析：一个研究者也可使用两个以上的自变量，这是一个N向方差分析（N是你拥有的自变量的数量），也就是MANOVA检验。

    • 例如，能够同时按国家、性别、年龄组、种族等检查日冕病例的潜在差别
    • 方差分析会给你一个单变量的f值，而方差分析会给你一个多变量的f值

有复制与无复制

你可能常常听到关于方差分析的复制和不复制。让咱们了解这些是什么：

  1. 具备复制功能的双向ANOVA：两个小组和这些小组的成员所作的不仅是一件事情

    • 例如，假设还没有开发出针对COVID-19的疫苗，医生正在尝试两种不一样的治疗方法来治愈两组感染COVID-19的患者

  2. 双向ANOVA（无复制）：只有一个组而且对同一组进行双重测试时使用

    • 例如，假设已为COVID-19开发了一种疫苗，研究人员正在对一组志愿者进行疫苗接种以前和以后的测试，以查看其是否有效

过后检验

当咱们进行方差分析时，咱们试图肯定各组之间是否存在统计学上的显著差别。若是咱们发现存在差别，则须要检查组差别的位置。

基本上，过后检验告诉研究者哪些组彼此不一样。

此时，你能够运行过后检验，这是t检验，用于检验组之间的均值差别。能够进行多个比较测试来控制I型错误率，包括Bonferroni、Scheffe、Dunnet和Tukey测试。

如今，让咱们用一些真实的数据来理解每种类型的方差分析测试，并使用Python。

Python中的单向方差分析测试

我从一个正在进行的Kaggle竞赛中下载了这些数据:https://www.kaggle.com/sudala...

在此测试中，咱们将尝试分析区域或状态的密度与日冕例数之间的关系。所以，咱们将根据每一个州的人口密度来映射每一个州。

首先导入全部必需的库和数据：

import pandas as pd
import numpy as np

import scipy.stats as stats
import os
import random

import statsmodels.api as sm
import statsmodels.stats.multicomp

from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns

从目录加载数据：

StatewiseTestingDetails=pd.read_csv('./StatewiseTestingDetails.csv')
population_india_census2011=pd.read_csv('./population_india_census2011.csv')

StatewiseTestingDetails包含有关每一个州一天中阳性和阴性病例总数的信息。而human_india_census2011包含有关每一个州的密度的信息以及有关人口的其余相关信息。

population_india_census2011.head() 
StatewiseTestingDetails.head() #了解数据
StatewiseTestingDetails['Positive'].sort_values().head() #排序

从上面的代码片断中，咱们能够看到有几个州在一天内有0个日冕案例或没有日冕案例。因此让咱们看看这样的州：

StatewiseTestingDetails['State'][StatewiseTestingDetails['Positive']==1].unique()

咱们看到，Nagaland和Sikkim在一天内也没有日冕病例。另外一方面，Arunachal和Mizoram一天只有一个日冕病例。

估算缺失值：咱们注意到“Positive”列中有许多缺失值。所以，让咱们用每一个州的Positive中值来估算这些缺失的值：

stateMedianData=StatewiseTestingDetails.groupby('State')[['Positive']].median().reset_index().rename(columns={'Positive':'Median'})
stateMedianData.head()

for index,row in StatewiseTestingDetails.iterrows():

    if pd.isnull(row['Positive']):

        StatewiseTestingDetails['Positive'][index]=int(stateMedianData['Median'][stateMedianData['State']==row['State']])

#合并StatewiseTestingDetails & population_india_census2011
data=pd.merge(StatewiseTestingDetails,population_india_census2011,on='State')

如今咱们能够编写一个函数，根据每一个州的密度建立一个密度组bucket，其中Dense1<Dense2<Dense3<Dense4：

def densityCheck(data):
    data['density_Group']=0
    for index,row in data.iterrows():
        status=None
        i=row['Density'].split('/')[0]
        try:
            if (',' in i):
                i=int(i.split(',')[0]+i.split(',')[1])
            elif ('.' in i):
                i=round(float(i))
            else:
                i=int(i)
        except ValueError as err:
            pass
        try:
            if (0<i<=300):
                status='Dense1'
            elif (300<i<=600):
                status='Dense2'
            elif (600<i<=900):
                status='Dense3'
            else:
                status='Dense4'
        except ValueError as err:
            pass
        data['density_Group'].iloc[index]=status
    return data

如今，用密度组映射每一个州。咱们能够导出这些数据，以便之后在双向方差分析测试中使用：

data=densityCheck(data)
#导出
stateDensity=data[['State','density_Group']].drop_duplicates().sort_values(by='State')

让咱们对能够用于方差分析测试的数据集进行从新排列：

df=pd.DataFrame({'Dense1':data[data['density_Group']=='Dense1']['Positive'],
                 'Dense2':data[data['density_Group']=='Dense2']['Positive'],
                 'Dense3':data[data['density_Group']=='Dense3']['Positive'],
                 'Dense4':data[data['density_Group']=='Dense4']['Positive']})

咱们的ANOVA检验假设之一是应随机选择样本，而且样本应接近高斯分布。所以，让咱们从每一个因子或水平中选择10个随机样本：

np.random.seed(1234)
dataNew=pd.DataFrame({'Dense1':random.sample(list(data['Positive'][data['density_Group']=='Dense1']), 10),
'Dense2':random.sample(list(data['Positive'][data['density_Group']=='Dense1']), 10),
'Dense3':random.sample(list(data['Positive'][data['density_Group']=='Dense1']), 10),
'Dense4':random.sample(list(data['Positive'][data['density_Group']=='Dense1']), 10)})

让咱们绘制日冕案例数量的密度分布图，以检查它们在不一样密度组中的分布：

咱们能够清楚地看到数据不遵循高斯分布。

有不一样的数据转换方法可使数据接近高斯分布。咱们进行Box-Cox变换：

dataNew['Dense1'],fitted_lambda = stats.boxcox(dataNew['Dense1'])
dataNew['Dense2'],fitted_lambda = stats.boxcox(dataNew['Dense2'])
dataNew['Dense3'],fitted_lambda = stats.boxcox(dataNew['Dense3'])
dataNew['Dense4'],fitted_lambda = stats.boxcox(dataNew['Dense4'])

如今让咱们再次绘制它们的分布图来检查：

方法1：使用statsmodels模块进行单向方差分析

Python中有两种方法能够执行ANOVA测试。一个是stats.f_oneway()方法：

F, p = stats.f_oneway(dataNew['Dense1'],dataNew['Dense2'],dataNew['Dense3'],dataNew['Dense4'])
## 看看总体模型是否重要
print('F-Statistic=%.3f, p=%.3f' % (F, p))

咱们发现p值<0.05。所以，咱们能够拒绝零假设——不一样密度组之间没有差别。

方法2：用OLS模型进行单因素方差分析

正如咱们在回归中所知道的，咱们能够对每一个输入变量进行回归，并检查其对目标变量的影响。因此，咱们将遵循一样的方法，咱们在线性回归中遵循的方法。

model = ols('Count ~ C(density_Group)', newDf).fit()
model.summary()

## 看看总体模型是否重要
print(f"Overall model F({model.df_model: .0f},{model.df_resid: .0f}) = {model.fvalue: .3f}, p = {model.f_pvalue: .4f}")

#建立方差分析表
res = sm.stats.anova_lm(model, typ= 2)
res

从以上输出结果能够看出，p值小于0.05。所以，咱们能够拒毫不同密度组之间没有差别的零假设。

F-statistic= 5.817，p-value= 0.002，这代表density_Group对日冕阳性病例有整体显着影响。可是，咱们尚不知道desnity_groups之间的区别在哪里。所以，基于p值，咱们能够拒绝H0；就面积密度和日冕例数而言，没有显着差别。

过后比较检验

当咱们进行方差分析时，咱们试图肯定各组之间是否存在统计学上的显着差别。那么，若是咱们发现统计学意义呢？

若是发现存在差别，则须要检查组差别的位置。所以，咱们将使用Tukey HSD测试来肯定差别所在：

mc = statsmodels.stats.multicomp.MultiComparison(newDf['Count'],newDf['density_Group'])
mc_results = mc.tukeyhsd()
print(mc_results)

Tuckey HSD测试清楚地代表，Group1 – Group3，Group1 – Group4，Group2 – Group3和Group3 – Group4之间存在显着差别。

这代表，除上述两组外，全部其余日冕病例数的成对比较均拒绝零假设，且无统计学显著性差别。

假设检验/模型诊断

正态分布假设检验

当使用线性回归和方差分析模型时，假设与残差有关，而不是变量自己。

方法1:Shapiro-Wilk试验：
### 正态性假设检查
w, pvalue = stats.shapiro(model.resid)
print(w, pvalue)

从上面的代码片断中，咱们看到全部密度组的p值都大于0.05。所以，咱们能够得出结论，它们遵循高斯分布。

方法2:Q-Q图试验：

咱们可使用Q-Q图来检验这个假设：

res = model.resid
fig = sm.qqplot(res, line='s')
plt.show()

从上图中，咱们看到全部数据点都靠近45度线，所以咱们能够得出结论，它遵循正态分布。

方差假设检验的同质性检查

应针对分类变量的每一个级别检查方差假设的同质性。咱们可使用Levene检验来检验组之间的均等方差。

w, pvalue = stats.bartlett(newDf['Count'][newDf['density_Group']=='Dense1'], newDf['Count'][newDf['density_Group']=='Dense2']
                           , newDf['Count'][newDf['density_Group']=='Dense3'], newDf['Count'][newDf['density_Group']=='Dense4'])
print(w, pvalue)

## Levene 方差测试
stats.levene(dataNew['Dense1'],dataNew['Dense2'],dataNew['Dense3'],dataNew['Dense4'])

咱们发现全部密度组的p值都大于0.05。所以，咱们能够得出结论，各组具备相等的方差。

Python中的双向方差分析测试

一样，使用相同的数据集，咱们将试图了解一个地区或州的密度、人口年龄和日冕病例数量之间是否存在显著关系。所以，咱们将根据居住在其中的人口密度绘制每一个州的地图。

让咱们导入数据并检查是否存在任何数据歧义：

individualDetails=pd.read_csv('./individualDetails.csv',parse_dates=[2])
stateDensity=pd.read_csv('./stateDensity.csv')

从上面的代码片断中，咱们能够看到没有感染婴儿的记录。接下来，检查数据中是否缺乏值：

individualDetails.isna().sum()

print('Percentage of missing values in age & gender columns respectively :', \
      (individualDetails['age'].isna().sum()/individualDetails.shape[0])*100,'%',\
      (individualDetails['gender'].isna().sum()/individualDetails.shape[0])*100,'%')

咱们发如今年龄和性别栏中分别有超过91%和80%的条目丢失。因此咱们须要设计一种方法来估算它们。

在这里，我将以各州的性别中位数和各州的性别中位数估算年龄。所以，我将计算中位数和众数：

ageMedianPerState=individualDetails[~individualDetails['age'].isna()]
ageMedianPerState['age']=ageMedianPerState['age'].astype(str).astype(int)
ageMedianPerState=ageMedianPerState.groupby('State')[['age']].median().reset_index()
ageMedianPerState['age']=ageMedianPerState['age'].apply(lambda x:math.ceil(x))

#经过COVID-19查找每一个州的最常感染的性别
genderModePerState=individualDetails.groupby(['State'])['gender'].agg(pd.Series.mode).to_frame().reset_index()

#这没有得到有关整体性别的信息性别
genderModePerState=genderModePerState[genderModePerState['State']!='Arunachal Pradesh']
#如今在年龄和性别栏填充丢失的值
for index,row in individualDetails.iterrows():
    if row['State']=='Arunachal Pradesh':    
        individualDetails.drop([index],inplace=True)
        continue
    if pd.isnull(row['age']):
        individualDetails['age'][index]=list(ageMedianPerState['age'][ageMedianPerState['State']=='West Bengal'].values)[0]
    if pd.isnull(row['gender']):
        if len(genderModePerState['gender'][genderModePerState['State']==row['State']].values)>0:
            individualDetails['gender'][index]=(genderModePerState['gender'][genderModePerState['State']==row['State']].values[0])

如今，让咱们合并individualDetails和stateDensity数据帧，为咱们建立一个总体数据集：

data = pd.merge(individualDetails,stateDensity,on='State',how='left').reset_index(drop=True)

如今咱们能够建立年龄组桶：

data.dropna(subset=['density_Group'],inplace=True)
data.reset_index(drop=True,inplace=True)

合并数据以得到一个数据集，其中每一个人都映射了他们的年龄组和各自的州密度组：

patient_Count=data.groupby(['diagnosed_date','density_Group'])[['diagnosed_date']].count().\

                        rename(columns={'diagnosed_date':'Count'}).reset_index()

data=pd.merge(data,patient_Count,on=['diagnosed_date','density_Group'],how='inner')

newData=data.groupby(['density_Group','age_Group'])['Count'].apply(list).reset_index()

newData.head()

np.random.seed(1234)

AnovaData=pd.DataFrame(columns=['density_Group','age_Group','Count'])

for index,row in newData.iterrows():
    count=17
    tempDf=pd.DataFrame(index=range(0,count),columns=['density_Group','age_Group','Count'])

    tempDf['age_Group']=newData['age_Group'][index]

    tempDf['density_Group']=newData['density_Group'][index]

    tempDf['Count']=random.sample(list(newData['Count'][index]),count)

    AnovaData=pd.concat([AnovaData,tempDf],axis=0)

检查数据中Count列的分布，并使用箱线图方法检查数据中是否存在异常值：

plt.hist(AnovaData['Count'])
plt.show() sns.kdeplot(AnovaData['Count'],cumulative=False,bw=2)

咱们发如今咱们的数据中有许多异常值。甚至计数变量的分布也不是高斯分布。所以，咱们将使用Box-Cox变换方法来处理这种状况：

sns.boxplot(x='age_Group', y='Count',
                 data=AnovaData,
                 palette="colorblind")

AnovaData['Count']=AnovaData['Count'].astype(int)

## 数据转换
AnovaData['newCount'],fitted_lambda = stats.boxcox(AnovaData['Count'])

import matplotlib.pyplot as plt
sns.kdeplot(AnovaData['newCount'],cumulative=False,bw=2)

如今让咱们使用OLS模型来检验咱们的假设：

## 拟合OlS模型-方法1
model2 = ols('newCount ~ C(age_Group)+ C(density_Group)', AnovaData).fit()

print(f"Overall model F({model2.df_model: },{model2.df_resid: }) = {model2.fvalue: }, p = {model2.f_pvalue: }")
model2.summary()

## 建立方差分析表
res2 = sm.stats.anova_lm(model2, typ=2)
res2

#检验残差的正态分布
res = model2.resid
fig = sm.qqplot(res, line='s')
plt.show()

从上面的Q-Q图，咱们能够看到残差几乎是正态分布的（尽管在最末端的点能够被贴现）。所以，咱们能够得出结论，它知足方差分析检验的正态性假设。

#方法2-检查组之间的交互
formula = 'newCount ~ C(age_Group) *C(density_Group)'
model = ols(formula, AnovaData).fit()
model.summary()

aov_table = anova_lm(model, typ=2)
print(aov_table.round(4))

结果的解释

经过ANOVA分析得到的日冕病例数，年龄组和密度组以及相互做用的P值具备统计学意义（P <0.05）。咱们得出结论，density_Group的类型显着影响日冕的结果。

age_Group显着影响日冕病例的结果，age_Group和density_Group的相互做用也显着影响日冕病例的结果。

过后检验

最后，让咱们肯定哪些组在统计上是不一样的。咱们将使用Tuckey HSD方法：

mc = statsmodels.stats.multicomp.MultiComparison(AnovaData['newCount'],AnovaData['density_Group'])
mc_results = mc.tukeyhsd()
print(mc_results)

mc = statsmodels.stats.multicomp.MultiComparison(AnovaData['newCount'],AnovaData['age_Group'])
mc_results = mc.tukeyhsd()
print(mc_results)

从上面的Tuckey HSD测试结果中，咱们能够清楚地看到，密度组中的Group1 – Group3，Group1 – Group4与年龄组中的Young – Adult＆Young –old组之间也存在显着差别。

所以，Tukey HSD的上述结果代表，除上述组外，日冕病例数的全部其余成对比较均拒绝了原假设，而且代表没有统计学上的显着差别。

结尾

在病毒大流行时期，我试着用一个相关的案例来解释方差分析。你能够克隆个人Github存储库来下载所有代码和数据：https://github.com/Praveen76/...

原文连接：https://www.analyticsvidhya.c...

欢迎关注磐创AI博客站：
http://panchuang.net/

sklearn机器学习中文官方文档：
http://sklearn123.com/

欢迎关注磐创博客资源汇总站：
http://docs.panchuang.net/

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程