作者|GUEST BLOG
編譯|VK
來源|Analytics Vidhya
介紹
“事實是每個人都相信的簡單陳述,也就是事實是沒有錯的,除非它被人發現了錯誤,假設有一個沒人愿意相信的建議,那么它要直到被發現有效的時候才能成為事實,” –愛德華·泰勒
我們正在應對一場空前規模的流行病,全世界的研究人員都在瘋狂地試圖開發一種疫苗或COVID-19的治療方法,而醫生們正試圖阻止這種流行病席卷整個世界,
我最近有了一個想法,把我的統計知識應用到這些大量COVID資料中,

考慮這樣一個場景:醫生有四種醫療方法來治療病人,一旦我們有了測驗結果,用最少時間治愈病人的治療會是最好的方法,
但如果這些病人中的一些已經部分治愈,或者其他藥物已經在治療他們呢?
為了作出一個有信心和可靠的決定,我們需要證據來支持我們的做法,這就是方差分析的概念發揮作用的地方,
在本文中,我將向你介紹方差分析測驗及其用于做出更好決策的不同型別,我將在Python中演示每種型別的ANOVA(方差分析)測驗,以可視化它們并處理COVID-19資料,
注意:你必須了解統計學的基本知識才能理解這個主題,最好了解t檢驗和假設檢驗,
什么是方差分析測驗(ANOVA)
方差分析,或稱方差分析,可以看作是兩組以上的t檢驗的推廣,獨立t檢驗用于比較兩組之間的條件平均值,當我們想比較兩組以上患者的病情平均值時,使用方差分析,
方差分析測驗模型中某個地方的平均值是否存在差異(測驗是否存在整體效應),但它不能告訴我們差異在哪里(如果存在),為了找出兩組之間的區別,我們必須進行事后檢驗,
要執行任何測驗,我們首先需要定義原假設和替代假設:
- 零假設–各組之間無顯著差異
- 替代假設–各組之間存在顯著差異
基本上,方差分析是通過比較兩種型別的變化來完成的,即樣本均值之間的變化,以及每個樣本內部的變化,以下公式表示單向Anova測驗統計資料,
ANOVA公式的結果,即F統計量(也稱為F比率),允許對多組資料進行分析,以確定樣本之間和樣本內部的可變性,
單向ANOVA的公式可以這樣寫:


當我們繪制ANOVA表時,上面的所有組成部分都可以如下所示:

一般來說,如果與F相關聯的p值小于0.05,則將拒絕原假設并支持替代假設,如果原假設被拒絕,我們可以得出結論,所有組的均值不相等,
注:如果被測組之間不存在真正的差異,也就是所謂的零假設,那么方差分析的F比統計結果將接近1,
ANOVA檢驗的假設
在進行方差分析之前,我們需要做一些假設:
- 從因子水平定義的總體中獨立且隨機地獲得觀察結果
- 每個因子水平的資料均呈正態分布
- 案例獨立性:樣本案例應相互獨立
- 方差的同質性:同質性是指各組之間的方差應近似相等
方差同質性的假設可以用Levene檢驗或Brown-Forsythe檢驗來檢驗,分數分布的正態性可以用直方圖、偏度和峰度值來檢驗,也可以用Shapiro-Wilk或Kolmogorov-Smirnov或Q-Q圖來檢驗,獨立性的假設可以根據研究設計來確定,
值得注意的是,方差分析對于假設獨立性的違規行為并不強大,這就是說,即使你違反了同質性或正態性的假設,你也可以進行測驗并基本相信結果,
但是,如果違反了獨立性假設,方差分析的結果是無效的,一般來說,在違反同質性的情況下,如果具有相同大小的組,則分析被認為是可靠的,對于違反正態性的情況,如果樣本量較大,繼續進行方差分析通常是可以的,
方差分析檢驗型別
-
單向方差分析:單向方差分析只有一個自變數
- 例如,可以按國家/地區評估日冕案例的差異,并且一個國家可以將2個,20個或更多不同的類別進行比較
-
雙向方差分析:雙向方差分析(也稱為因子方差分析)是指使用兩個獨立變數的方差分析
- 擴展上面的示例,雙向方差分析可以按年齡組(獨立變數1)和性別(獨立變數2)檢查日冕病例(因變數)的差異,雙向方差分析可用于檢查兩個獨立變數之間的相互作用,相互作用表明,自變數的所有類別之間的差異不是統一的
- 例如,老年組總體上可能比青年組具有更高的日冕病例,但是與歐洲國家相比,亞洲國家的差異可能更大(或更小)
-
N向方差分析:一個研究者也可以使用兩個以上的自變數,這是一個N向方差分析(N是你擁有的自變數的數量),也就是MANOVA檢驗,
- 例如,可以同時按國家、性別、年齡組、種族等檢查日冕病例的潛在差異
- 方差分析會給你一個單變數的f值,而方差分析會給你一個多變數的f值
有復制與無復制
你可能經常聽到關于方差分析的復制和不復制,讓我們了解這些是什么:
-
具有復制功能的雙向ANOVA:兩個小組和這些小組的成員所做的不只是一件事情
- 例如,假設尚未開發出針對COVID-19的疫苗,醫生正在嘗試兩種不同的治療方法來治愈兩組感染COVID-19的患者
-
雙向ANOVA(無復制):只有一個組并且對同一組進行雙重測驗時使用
- 例如,假設已為COVID-19開發了一種疫苗,研究人員正在對一組志愿者進行疫苗接種之前和之后的測驗,以查看其是否有效
事后檢驗
當我們進行方差分析時,我們試圖確定各組之間是否存在統計學上的顯著差異,如果我們發現存在差異,則需要檢查組差異的位置,
基本上,事后檢驗告訴研究者哪些組彼此不同,
此時,你可以運行事后檢驗,這是t檢驗,用于檢驗組之間的均值差異,可以進行多個比較測驗來控制I型錯誤率,包括Bonferroni、Scheffe、Dunnet和Tukey測驗,
現在,讓我們用一些真實的資料來理解每種型別的方差分析測驗,并使用Python,
Python中的單向方差分析測驗
我從一個正在進行的Kaggle競賽中下載了這些資料:https://www.kaggle.com/sudalairajkumar/covid19-in-india
在此測驗中,我們將嘗試分析區域或狀態的密度與日冕例數之間的關系,因此,我們將根據每個州的人口密度來映射每個州,
首先匯入所有必需的庫和資料:
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=https://www.cnblogs.com/panchuangai/p/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=https://www.cnblogs.com/panchuangai/p/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=https://www.cnblogs.com/panchuangai/p/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= https://www.cnblogs.com/panchuangai/p/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 = https://www.cnblogs.com/panchuangai/p/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 = https://www.cnblogs.com/panchuangai/p/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 = https://www.cnblogs.com/panchuangai/p/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=https://www.cnblogs.com/panchuangai/p/pd.merge(data,patient_Count,on=['diagnosed_date','density_Group'],how='inner')
newData=https://www.cnblogs.com/panchuangai/p/data.groupby(['density_Group','age_Group'])['Count'].apply(list).reset_index()
newData.head()
np.random.seed(1234)
AnovaData=https://www.cnblogs.com/panchuangai/p/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=https://www.cnblogs.com/panchuangai/p/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=https://www.cnblogs.com/panchuangai/p/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/ANOVA-Test-COVID-19
原文鏈接:https://www.analyticsvidhya.com/blog/2020/06/introduction-anova-statistics-data-science-covid-python/
歡迎關注磐創AI博客站:
http://panchuang.net/
sklearn機器學習中文官方檔案:
http://sklearn123.com/
歡迎關注磐創博客資源匯總站:
http://docs.panchuang.net/
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/14686.html
標籤:其他
上一篇:Tensorflow 2.0 DLL load failed: 找不到指定的模塊
下一篇:深入淺出Transformer
