目錄

前言

為什么不能兩兩比較？方差分析（ANOVA）原理實例講解3.1提出問題3.2畫圖觀察3.3計算各誤差平方和3.4計算F檢驗值3.5R語言代碼判定系數事后檢驗參考資料后記前言

我們知道，在比較兩個分組之間有沒有差異時，我們會首選Ttest進行分析。如果樣本量太小或者數據分布不滿足正態性時，我們會選擇[Wilcoxon檢驗]Wilcoxon檢驗-簡書(jianshu.com)。

但是，在我們課題中，我們的實驗組可能不止2組，例如：用A藥組+用B藥組+用C藥組+用D藥組+……

在這種情況下，我們該怎么辦呢？

1.為什么不能兩兩比較？

最簡單來說，我們可能會想著把所有的組分成兩兩比較的組。例如：對于A組+B組+C組+D組來說，我們可能會用Ttest去分別比較：

A+B

A+C

A+D

B+C

B+D

C+D但是這里會存在一些問題：

檢驗過程繁瑣：比較次數多，這里因為只有4組，所以一共需要比較C42=6次Ttest進行檢驗。當分組更多時，檢驗次數則呈現指數化的增長。

檢驗的靈敏度降低：每次比較只用了部分數據，沒有考慮整體數據變化，使得計算結果的準確性降低，從而降低了檢驗的靈敏度。

結果可靠性低：在我們認同的a=0.05的前提下，也就是每次檢驗正確的可能性是0.95。那么6次檢驗同時正常的概率=(0.95)6=0.735。那么這個時候檢驗的a=1-0.735=0.265，遠遠大于我們普遍接受的0.05

所以我們需要一種檢驗來同時比較多組均值之間的關系，這就引出了我們的方差分析檢驗

2.方差分析（ANOVA）原理1.方差分析的簡介

它的基本思想是將測量數據的總變異（即總方差）按照變異來源分為處理（組間）效應和誤差（組內）效應，并作出其數量估計，從而確定實驗處理對研究結果影響力的大小。

按因素劃分，可分為單因素方差分析、二因素方差分析和多因素方差分析。

方差分析的步驟為：總平方和分解→總自由度分解→F檢驗。若F檢驗顯著，則可以進行多重比較，從而發現哪些處理具有兩兩間的差異。

2.方差分析的基本原理

在一次實驗中，可以得到一系列不同的觀測值。造成觀測值不同的原因可能是①由于處理因素不同引起的，即處理效應；也可能是②由于實驗過程中偶然性因素的干擾和測量誤差所致，即誤差效應。反應測量數據變異性的指標有多個，在方差分析中選用方差來度量資料的變異程度。要正確認識觀測值的變異是由于處理效應還是誤差效應引起的，我們可以分別計算出處理效應的方差以及誤差效應的方差，在一定顯著水平下進行比較，如果二者相差不大，說明實驗處理對觀測值的影響不大；如果差異較大則說明實驗處理對于觀測的影響較大。

全體樣本的總誤差平方和（總變異），也成為SST：[圖片上傳失敗每個組內的誤差平方和（組內變異），也稱為SSE：[圖片上傳失敗每個組間的誤差平方和（組間變異），也稱為SSA：[圖片上傳失敗他們的之間的數學關系：SST=SSE+SSA

i：i組，例子中的i可以理解為a,b,c,d

nj：第j組的總人數

j：每個組中的單個樣本

：全體樣本數據的總均值

：每組樣本的均值

：第i組中的第j個樣本

image.png

我們把變異分為組內變異和組間變異兩部分：

image.png

仔細觀察下面的例子，左圖各組樣本均數不相同，而右圖則樣本均數一致。

那我們再看下其兩類組間變異：左圖中組間差異很大，右圖中組間差異很小。

我們再看下其兩類組內變異：，左圖中組內差異較小，而右圖中組內差異較大。

image.png

如果組間差異（B）遠大于組內差異（A），就意味著各組樣本均數不一致呢？

方差分析就是基于這樣的思路：

以組內差異（A）為參考基準，考察組間差異（B）的大小。

如果組間差異（B）遠大于組內差異（A），則認為組間存在區別。

而組內差異，我們認為是因為（完全）隨機而產生的。以這樣一個完全隨機的尺度作為標桿，也甚是巧妙。

我們也引出了F值，即B比A的值。

image.png

基于F分布，就很容易看出，當組間差異越大，那么F=B/A值也越大，橫坐標越向右，越容易進入我們拒絕原假設（H0，各組均數相同）的區域。

也就是說，當組間差異越大，p值越小，越有理由拒絕原假設H0，也就是兩組之間的差異越有顯著性。

3.實例講解3.1提出問題

假設我們收集了不同藥物降低體重的數據：

image.png

那么這4種藥物之間的效果有沒有差異呢？

3.2畫圖觀察

我們首先畫個散點圖看看，具體畫圖代碼如下：

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)

結果如下：

image.png

從圖中我們可以大致看出這4組藥物的效果之間存在一定差異。

3.3計算各誤差平方和

SST的計算：直接計算即可，計算得到4164.609

SSA的計算：直接計算即可，計算得到1456.609

SSE的計算：直接計算即可，計算得到2708

這些指標的計算代碼如下：

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#定義函數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

這里再次印證了前面的結論：SST=SSE+SSA

3.4計算F檢驗值

1.因為誤差平方和會受到樣本數的影響，所以為了消除樣本數的影響，我們將各個誤差平方和除以（樣本數-1），其實也就是方差

2.綜合前面所有數據，得到的結果值統計如下：

image.png

(*)：此處的方差不可通過MSA+MSE計算而來

3.如果不同組藥物之間的效果一致，或者相差不大，那么MSA和MSE之間的比值與1接近。所以我們需要對F-test=MSA/MSE進行檢驗，得到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

可以看到，這個結果和我們手動計算的結果是一致的！F=3.407，對應的P=0.0388

4.判定系數

1.判定系數（CoefficientofDetermination，記為R2）在統計學中用于度量因變量的變異中可由自變量解釋部分所占的比例，以此來判斷統計模型的解釋力。

2.決定系數是相關系數的二次冪，例如在這個例子中，我們只有一個自變量——不同的治療方案。通過計算不同方案的決定系數，我們可以得知不同治療方案和效果的相關系數R。

3.決定系數意義：擬合優度越大，自變量對因變量的解釋程度越高，自變量引起的變動占總變動的百分比高。觀察點在回歸直線附近越密集。

4.在方差分析中，我們可以通過計算全部變異（SST）中，由分組不同所解釋變異（SSA）的占比（SSA/SST）來計算判定系數：

image.png

所以可以計算得知相關系數R：

image.png

上述數據說明了：

R2=0.3498：不同的藥物解釋了數據34.98%的變異。

R=0.5914：不同的藥物選擇與療效之間的相關性為0.5914。

5.事后檢驗

1.有了方差分析的結果，我們只能說明j個總體均值不全相等。若想進一步了解哪些兩個總體均值不等，需進行多個樣本均值間的兩兩比較或稱多重比較（multiplecomparison），也稱為posthoc（事后）檢驗。

2.事后檢驗的***非常多，這里不展開說明，僅列出常見的***：

SNK-q檢驗

LSD-t檢驗：同過檢驗兩個總體均值是否相等的T檢驗，其中總體方差用MSE來代替而得到的。

image.png

Dunnett檢驗

TukeyHSD檢驗

Duncan檢驗

Scheffe檢驗1.下面我用TukeyHSD檢驗對我們的結果進行檢驗，并告訴大家如何解讀自己的結果。這里只展示代碼，具體原理日后有機會再補充。

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

可以看到，在事后檢驗中，僅僅只有d-c之間存在顯著差異，而其他檢驗之間沒有顯著差異。

我們可以把結果進行繪制圖形：

1plot(TukeyHSD(fit))