【原】十步搞定單因素方差分析

健明 2021-07-14

展開全文

看過(guò)我TCGA腫瘤數(shù)據(jù)庫(kù)知識(shí)圖譜的小伙伴都只是如何在任意癌癥查詢指定感興趣基因的表達(dá)量，并且對(duì)樣本進(jìn)行分組比較，網(wǎng)站是：https:///heatmap/
悄咪咪的上線了TCGA知識(shí)圖譜視頻教程（B站和YouTube直達(dá)）

image-20190506132847002

根據(jù)視頻教程拿到數(shù)據(jù)，很容易可視化如下：

image-20190506132955029

可以得出結(jié)論，我們感興趣的基因（這里是CUL5）在乳腺癌的正常組織及癌癥組織（原位和轉(zhuǎn)移）表達(dá)量，使用單因素方差分析，得到了統(tǒng)計(jì)學(xué)顯著的結(jié)果。

定義

單因素方差分析是兩個(gè)樣本平均數(shù)比較的引伸，它是用來(lái)檢驗(yàn)多個(gè)平均數(shù)之間的差異，從而確定因素對(duì)試驗(yàn)結(jié)果有無(wú)顯著性影響的一種統(tǒng)計(jì)方法。

因素：影響研究對(duì)象的某一指標(biāo)、變量。
水平：因素變化的各種狀態(tài)或因素變化所分的等級(jí)或組別。
單因素試驗(yàn)：考慮的因素只有一個(gè)的試驗(yàn)叫單因素試驗(yàn)。

了解數(shù)據(jù)

數(shù)據(jù)文件可以下載，然后讀入R里面進(jìn)行可視化，代碼如下：

rm(list = ls())
options(stringsAsFactors = F)
library(ggpubr)
df=read.table('CUL5-BRCA-type.tsv',header = T,sep = '\t')
colnames(df)

p = ggboxplot(df, "sample_type", "CUL5",
          color = "sample_type", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
          add = "jitter", shape = "sample_type")
p

這個(gè)文件 CUL5-BRCA-type.tsv 如果你沒(méi)有看我的TCGA腫瘤數(shù)據(jù)庫(kù)知識(shí)圖譜可能不知道如何下載，可以發(fā)郵件給我找我申請(qǐng)這個(gè)測(cè)試數(shù)據(jù) ( 郵箱： jmzeng1314@163.com )

可以看到，比網(wǎng)頁(yè)工具出圖要好看：

image-20190506133842005

但是，這個(gè)時(shí)候還沒(méi)有進(jìn)行統(tǒng)計(jì)分析，可以添加的統(tǒng)計(jì)學(xué)檢驗(yàn)包括：

image-20190506133655190

代碼也很簡(jiǎn)單：

p+stat_compare_means(method = "anova", label.y = 10)+      # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.")                  # Pairwise comparison against all

image-20190506134047019

可以看到，跟網(wǎng)頁(yè)工具結(jié)果一模一樣，而且出圖更漂亮，下面我們就手把手帶領(lǐng)大家完成這個(gè)分析，把這個(gè)ggpubr一步就完成的工作拆解開來(lái)。

第1步：計(jì)算各組內(nèi)樣本均值

lapply(split(df,df$sample_type),function(x) mean(x$CUL5))

第2步：計(jì)算所有樣本均值

mean(df$CUL5 )

第3步：計(jì)算各組內(nèi)部誤差平方和

tmp=lapply(split(df,df$sample_type),function(x) sum((x$CUL5-mean(x$CUL5))^2) )
sse = sum(unlist(tmp))

第4步：計(jì)算各組間誤差平方和

tmp=lapply(split(df,df$sample_type),function(x) nrow(x)*(mean(x$CUL5) - mean(df$CUL5 ))^2 )
ssb = sum(unlist(tmp))

第5步：計(jì)算各組內(nèi)部均方誤

mse=sse/(nrow(df)-length(unique(df$sample_type)))

第6步：計(jì)算組間均方誤

msd=ssb/length(unique(df$sample_type))-1

第7步：計(jì)算F比率

f= msb/mse
f

第8步：查找F臨界值

df1=(length(unique(df$sample_type))-1)
df2=(nrow(df)-length(unique(df$sample_type)))
qf(0.05,2,1215)

差別可知這里的F值是0.05，遠(yuǎn)小于我們真實(shí)情況，所以非常顯著了。

第9步：判斷是否顯著

1-pf(f,2,1215)

現(xiàn)在我們已經(jīng)知道了，在選定的顯著水平為0.05時(shí)候，這個(gè)F統(tǒng)計(jì)是顯著的，但是仍然是不知道哪組之間不一樣，所以可以選擇tukey檢驗(yàn)

第10步：進(jìn)行tukey檢驗(yàn),多重比較

J·W·圖凱（Tukey）于1953年提出一種能將所有各對(duì)平均值同時(shí)比較的方法，這種方法現(xiàn)在已被廣泛采用，一般稱之為“HSD檢驗(yàn)法”，或稱“W法”。 Tukey (John Wilder Tukey) for multiple comparisons
主要應(yīng)用于3組或以上的多重比較。比如說(shuō)一共有4組數(shù)據(jù)，兩兩比較產(chǎn)生6個(gè)統(tǒng)計(jì)值，Tukey test用于生成一個(gè)critical value來(lái)控制總體誤差（Familywise error rate，F(xiàn)ER）;與Tukey test相類似的是Dunnett test，它是控制多對(duì)一比較（即3組同時(shí)和一個(gè)參照組比較）的FER。

這個(gè)多重比較算法還蠻多的，參考：https://zhuanlan.zhihu.com/p/44880434 這里超綱了，我就不具體介紹了！