方差分析介紹(結合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教程