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

健明 2021-07-14

展開全文

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

image-20190506132847002

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

image-20190506132955029

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

定義

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

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

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

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

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

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

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

image-20190506133842005

但是，這個時候還沒有進行統(tǒng)計分析，可以添加的統(tǒng)計學(xué)檢驗包括：

image-20190506133655190

代碼也很簡單：

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)頁工具結(jié)果一模一樣，而且出圖更漂亮，下面我們就手把手帶領(lǐng)大家完成這個分析，把這個ggpubr一步就完成的工作拆解開來。

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

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

第2步：計算所有樣本均值

mean(df$CUL5 )

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

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

第4步：計算各組間誤差平方和

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

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

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

第6步：計算組間均方誤

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

第7步：計算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，遠小于我們真實情況，所以非常顯著了。

第9步：判斷是否顯著

1-pf(f,2,1215)

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

第10步：進行tukey檢驗,多重比較

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

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