前言
為什么不能兩兩比較？方差分析（ANOVA）原理實(shí)例講解3.1提出問題3.2畫圖觀察3.3計(jì)算各誤差平方和3.4計(jì)算F檢驗(yàn)值3.5R語言代碼判定系數(shù)事后檢驗(yàn)參考資料后記前言我們知道，在比較兩個分組之間有沒有差異時，我們會首選Ttest進(jìn)行分析。如果樣本量太小或者數(shù)據(jù)分布不滿足正態(tài)性時，我們會選擇[Wilcoxon檢驗(yàn)]Wilcoxon檢驗(yàn)-簡書(jianshu.com)。
但是，在我們課題中，我們的實(shí)驗(yàn)組可能不止2組，例如：用A藥組+用B藥組+用C藥組+用D藥組+……
在這種情況下，我們該怎么辦呢？
1.為什么不能兩兩比較？最簡單來說，我們可能會想著把所有的組分成兩兩比較的組。例如：對于A組+B組+C組+D組來說，我們可能會用Ttest去分別比較：
A+B
A+C
A+D
B+C
B+D
C+D但是這里會存在一些問題：
檢驗(yàn)過程繁瑣：比較次數(shù)多，這里因?yàn)橹挥?組，所以一共需要比較C42=6次Ttest進(jìn)行檢驗(yàn)。當(dāng)分組更多時，檢驗(yàn)次數(shù)則呈現(xiàn)指數(shù)化的增長。
檢驗(yàn)的靈敏度降低：每次比較只用了部分?jǐn)?shù)據(jù)，沒有考慮整體數(shù)據(jù)變化，使得計(jì)算結(jié)果的準(zhǔn)確性降低，從而降低了檢驗(yàn)的靈敏度。
結(jié)果可靠性低：在我們認(rèn)同的a=0.05的前提下，也就是每次檢驗(yàn)正確的可能性是0.95。那么6次檢驗(yàn)同時正常的概率=(0.95)6=0.735。那么這個時候檢驗(yàn)的a=1-0.735=0.265，遠(yuǎn)遠(yuǎn)大于我們普遍接受的0.05
所以我們需要一種檢驗(yàn)來同時比較多組均值之間的關(guān)系，這就引出了我們的方差分析檢驗(yàn)
2.方差分析（ANOVA）原理1.方差分析的簡介它的基本思想是將測量數(shù)據(jù)的總變異（即總方差）按照變異來源分為處理（組間）效應(yīng)和誤差（組內(nèi)）效應(yīng)，并作出其數(shù)量估計(jì)，從而確定實(shí)驗(yàn)處理對研究結(jié)果影響力的大小。
按因素劃分，可分為單因素方差分析、二因素方差分析和多因素方差分析。
方差分析的步驟為：總平方和分解→總自由度分解→F檢驗(yàn)。若F檢驗(yàn)顯著，則可以進(jìn)行多重比較，從而發(fā)現(xiàn)哪些處理具有兩兩間的差異。
2.方差分析的基本原理在一次實(shí)驗(yàn)中，可以得到一系列不同的觀測值。造成觀測值不同的原因可能是①由于處理因素不同引起的，即處理效應(yīng)；也可能是②由于實(shí)驗(yàn)過程中偶然性因素的干擾和測量誤差所致，即誤差效應(yīng)。反應(yīng)測量數(shù)據(jù)變異性的指標(biāo)有多個，在方差分析中選用方差來度量資料的變異程度。要正確認(rèn)識觀測值的變異是由于處理效應(yīng)還是誤差效應(yīng)引起的，我們可以分別計(jì)算出處理效應(yīng)的方差以及誤差效應(yīng)的方差，在一定顯著水平下進(jìn)行比較，如果二者相差不大，說明實(shí)驗(yàn)處理對觀測值的影響不大；如果差異較大則說明實(shí)驗(yàn)處理對于觀測的影響較大。
全體樣本的總誤差平方和（總變異），也成為SST：[圖片上傳失敗每個組內(nèi)的誤差平方和（組內(nèi)變異），也稱為SSE：[圖片上傳失敗每個組間的誤差平方和（組間變異），也稱為SSA：[圖片上傳失敗他們的之間的數(shù)學(xué)關(guān)系：SST=SSE+SSAi：i組，例子中的i可以理解為a,b,c,d
nj：第j組的總?cè)藬?shù)
j：每個組中的單個樣本
：全體樣本數(shù)據(jù)的總均值
：每組樣本的均值
：第i組中的第j個樣本
image.png
我們把變異分為組內(nèi)變異和組間變異兩部分：
image.png
仔細(xì)觀察下面的例子，左圖各組樣本均數(shù)不相同，而右圖則樣本均數(shù)一致。
那我們再看下其兩類組間變異：左圖中組間差異很大，右圖中組間差異很小。
我們再看下其兩類組內(nèi)變異：，左圖中組內(nèi)差異較小，而右圖中組內(nèi)差異較大。
image.png
如果組間差異（B）遠(yuǎn)大于組內(nèi)差異（A），就意味著各組樣本均數(shù)不一致呢？
方差分析就是基于這樣的思路：
以組內(nèi)差異（A）為參考基準(zhǔn)，考察組間差異（B）的大小。
如果組間差異（B）遠(yuǎn)大于組內(nèi)差異（A），則認(rèn)為組間存在區(qū)別。
而組內(nèi)差異，我們認(rèn)為是因?yàn)椋ㄍ耆╇S機(jī)而產(chǎn)生的。以這樣一個完全隨機(jī)的尺度作為標(biāo)桿，也甚是巧妙。
我們也引出了F值，即B比A的值。
image.png
基于F分布，就很容易看出，當(dāng)組間差異越大，那么F=B/A值也越大，橫坐標(biāo)越向右，越容易進(jìn)入我們拒絕原假設(shè)（H0，各組均數(shù)相同）的區(qū)域。
也就是說，當(dāng)組間差異越大，p值越小，越有理由拒絕原假設(shè)H0，也就是兩組之間的差異越有顯著性。
3.實(shí)例講解3.1提出問題假設(shè)我們收集了不同藥物降低體重的數(shù)據(jù)：
image.png
那么這4種藥物之間的效果有沒有差異呢？
3.2畫圖觀察我們首先畫個散點(diǎn)圖看看，具體畫圖代碼如下：
library(ggplot2)2dat<-data.frame(a=c(57,66,49,40,34,53,44),3b=c(68,39,29,45,56,51,NA),4c=c(31,49,21,34,40,NA,NA),5d=c(44,51,65,77,58,NA,NA))6gdat<-melt(dat,na.rm=T)7colnames(gdat)8gdat$variable<-as.character(gdat$variable)9mdat<-data.frame(x=c('a','b','c','d'),10y=sapply(dat,mean,na.rm=T))11ggplot()+12geom_point(gdat,mapping=aes(y=value,x=variable))+13geom_line(mdat,mapping=aes(y=y,x=x),color="red",group=1)結(jié)果如下：
image.png
從圖中我們可以大致看出這4組藥物的效果之間存在一定差異。
3.3計(jì)算各誤差平方和SST的計(jì)算：直接計(jì)算即可，計(jì)算得到4164.609
SSA的計(jì)算：直接計(jì)算即可，計(jì)算得到1456.609
SSE的計(jì)算：直接計(jì)算即可，計(jì)算得到2708
這些指標(biāo)的計(jì)算代碼如下：
1rm(list=ls())2options(stringsAsFactors=F)3library(reshape2)4dat<-data.frame(a=c(57,66,49,40,34,53,44),5b=c(68,39,29,45,56,51,NA),6c=c(31,49,21,34,40,NA,NA),7d=c(44,51,65,77,58,NA,NA))8gdat<-melt(dat,na.rm=T)9#定義函數(shù)diffsum10ds<-function(n,m){11sum((n-m)**2,na.rm=T)12}13#SST14ds(gdat$value,mean(gdat$value))#4164.60915#SSA16sum(((sapply(dat,mean,na.rm=T)-mean(gdat$value))**2)*table(gdat$variable))#1456.60917#SSE18sum(c(ds(dat$a,mean(dat$a,na.rm=T)),19ds(dat$b,mean(dat$b,na.rm=T)),20ds(dat$c,mean(dat$c,na.rm=T)),21ds(dat$d,mean(dat$d,na.rm=T))))#2708這里再次印證了前面的結(jié)論：SST=SSE+SSA
3.4計(jì)算F檢驗(yàn)值1.因?yàn)檎`差平方和會受到樣本數(shù)的影響，所以為了消除樣本數(shù)的影響，我們將各個誤差平方和除以（樣本數(shù)-1），其實(shí)也就是方差
2.綜合前面所有數(shù)據(jù)，得到的結(jié)果值統(tǒng)計(jì)如下：
image.png
(*)：此處的方差不可通過MSA+MSE計(jì)算而來3.如果不同組藥物之間的效果一致，或者相差不大，那么MSA和MSE之間的比值與1接近。所以我們需要對F-test=MSA/MSE進(jìn)行檢驗(yàn)，得到P值。
4.F-test值服從分子自由度為k-1、分母自由度為n-k的F分布，即：MSA/MSE=F-test~F(i-1,n-i)
5.查表得知P<0.05
3.5R語言代碼rm(list=ls())2options(stringsAsFactors=F)3library(reshape2)4dat<-data.frame(a=c(57,66,49,40,34,53,44),5b=c(68,39,29,45,56,51,NA),6c=c(31,49,21,34,40,NA,NA),7d=c(44,51,65,77,58,NA,NA))8gdat<-melt(dat,na.rm=T)9fit<-aov(value~variable,data=gdat)10summary(fit)11#DfSumSqMeanSqFvaluePr(>F)12#variable31457485.53.4070.0388*13#Residuals192708142.514#---15#Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1可以看到，這個結(jié)果和我們手動計(jì)算的結(jié)果是一致的！F=3.407，對應(yīng)的P=0.0388
4.判定系數(shù)1.判定系數(shù)（CoefficientofDetermination，記為R2）在統(tǒng)計(jì)學(xué)中用于度量因變量的變異中可由自變量解釋部分所占的比例，以此來判斷統(tǒng)計(jì)模型的解釋力。
2.決定系數(shù)是相關(guān)系數(shù)的二次冪，例如在這個例子中，我們只有一個自變量——不同的治療方案。通過計(jì)算不同方案的決定系數(shù)，我們可以得知不同治療方案和效果的相關(guān)系數(shù)R。
3.決定系數(shù)意義：擬合優(yōu)度越大，自變量對因變量的解釋程度越高，自變量引起的變動占總變動的百分比高。觀察點(diǎn)在回歸直線附近越密集。
4.在方差分析中，我們可以通過計(jì)算全部變異（SST）中，由分組不同所解釋變異（SSA）的占比（SSA/SST）來計(jì)算判定系數(shù)：
image.png
所以可以計(jì)算得知相關(guān)系數(shù)R：
image.png
上述數(shù)據(jù)說明了：
R2=0.3498：不同的藥物解釋了數(shù)據(jù)34.98%的變異。
R=0.5914：不同的藥物選擇與療效之間的相關(guān)性為0.5914。
5.事后檢驗(yàn)1.有了方差分析的結(jié)果，我們只能說明j個總體均值不全相等。若想進(jìn)一步了解哪些兩個總體均值不等，需進(jìn)行多個樣本均值間的兩兩比較或稱多重比較（multiplecomparison），也稱為posthoc（事后）檢驗(yàn)。
2.事后檢驗(yàn)的***非常多，這里不展開說明，僅列出常見的***：
SNK-q檢驗(yàn)
LSD-t檢驗(yàn)：同過檢驗(yàn)兩個總體均值是否相等的T檢驗(yàn)，其中總體方差用MSE來代替而得到的。
image.png
Dunnett檢驗(yàn)
TukeyHSD檢驗(yàn)
Duncan檢驗(yàn)
Scheffe檢驗(yàn)1.下面我用TukeyHSD檢驗(yàn)對我們的結(jié)果進(jìn)行檢驗(yàn)，并告訴大家如何解讀自己的結(jié)果。這里只展示代碼，具體原理日后有機(jī)會再補(bǔ)充。
1rm(list=ls())2options(stringsAsFactors=F)3library(reshape2)4dat<-data.frame(a=c(57,66,49,40,34,53,44),5b=c(68,39,29,45,56,51,NA),6c=c(31,49,21,34,40,NA,NA),7d=c(44,51,65,77,58,NA,NA))8gdat<-melt(dat,na.rm=T)9fit<-aov(value~variable,data=gdat)10TukeyHSD(fit)11#Tukeymultiplecomparisonsofmeans12#95%family-wiseconfidencelevel13#14#Fit:aov(formula=value~variable,data=gdat)15#16#$variable17#difflwruprpadj18#b-a-1-19.67609617.6760960.998739019#c-a-14-33.6560245.6560240.221814420#d-a10-9.65602429.6560240.496679021#c-b-13-33.3270707.3270700.304637822#d-b11-9.32707031.3270700.444784823#d-c242.76906845.2309320.0234078可以看到，在事后檢驗(yàn)中，僅僅只有d-c之間存在顯著差異，而其他檢驗(yàn)之間沒有顯著差異。
我們可以把結(jié)果進(jìn)行繪制圖形：
1plot(TukeyHSD(fit))
