R語言如何實現(xiàn)方差分析及其可視化

今天小編給大家分享一下R語言如何實現(xiàn)方差分析及其可視化的相關(guān)知識點，內(nèi)容詳細，邏輯清晰，相信大部分人都還太了解這方面的知識，所以分享這篇文章給大家參考一下，希望大家閱讀完這篇文章后有所收獲，下面我們一起來了解一下吧。

方差分析(Analysis of Variance，簡稱ANOVA)，又稱“變異數(shù)分析”，是R.A.Fisher發(fā)明的，用于兩個及兩個以上樣本均數(shù)差別的顯著性檢驗。 由于各種因素的影響，研究所得的數(shù)據(jù)呈現(xiàn)波動狀。 造成波動的原因可分成兩類，一是不可控的隨機因素，另一是研究中施加的對結(jié)果形成影響的可控因素。

方差分析分類：

方差分析常用的為單因素方差分析和雙因素方差分析，而多因素方差分析比較復雜用的不多，這里主要介紹前兩種方差分析。

單因素方差分析

單因素方差分析比較的是分類因子定義的兩個或多個組別中的因變量均值。與t檢驗區(qū)別：t檢驗適用于兩列數(shù)據(jù)的均值比較。單因素方差分析適用于兩列或更多列數(shù)據(jù)的均值比較。但對于兩列數(shù)據(jù)的均值比較，單因素方差分析=等方差假設(shè)的雙尾t檢驗。

以R中內(nèi)置數(shù)據(jù)PlantGrowth集為例。不同地塊（ ‘ctrl’, ‘trt1’, and ‘trt2’）種植的植物的產(chǎn)量（weight）數(shù)據(jù)；解決的問題如下
1.不同地塊產(chǎn)量是否存在差異
2.trt1與trt2是哪塊差異更大

##正態(tài)性檢驗

my_data <- PlantGrowth
shapiro.test(my_data$weight)

#結(jié)果表明，p=0.8915 >0.05，接受原假設(shè)，說明數(shù)據(jù)在多個水平下都是正態(tài)分布的。

##方差齊性檢驗
bartlett.test(weight~group,data=my_data)

#結(jié)果表明，p=0.2371 >0.05，接受原假設(shè)，說明數(shù)據(jù)在不同水平下是等方差的。

##單因素方差分析
fit<-aov(weight~group,data=my_data)
summary(fit)
#            Df Sum Sq Mean Sq F value Pr(>F)  
#group        2  3.766  1.8832   4.846 0.0159 *
#Residuals   27 10.492  0.3886                 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#結(jié)果表明，不同地塊見產(chǎn)量差異非常顯著。

#如果用非參數(shù)檢驗法 Kruskal-Wallis test
kruskal.test(weight ~ group, data = my_data)

##多重比較（多重T檢驗），哪塊差異更明顯
#R語言中TukeyHSD()函數(shù)提供了對各組均值差異的成對檢驗。
TukeyHSD(fit)


#如果用KW檢驗多重比較，
pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
                     p.adjust.method = "BH")

更好的方法+繪圖展示：

利用ggpubr包做方差分析，并可視化分析結(jié)果：

library("ggpubr")
compare_means(weight ~ group, data = PlantGrowth,method="anova")
compare_means(weight ~ group, data = PlantGrowth,method="kruskal.test")

my_comparisons <- list(c("ctrl","trt1"), c("ctrl", "trt2"),
                       c("trt1", "trt2"))

ggboxplot(my_data, x = "group", y = "weight", 
          color = "group", palette = "jco",
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment")+
          stat_compare_means(label.y = 11,method = "anova")+
          stat_compare_means(comparisons = my_comparisons,label.y = c(7,9,8))

或者小提琴圖

ggviolin(my_data, x = "group", y = "weight", 
          fill = "group", palette = "jco",
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment",add = "boxplot", add.params = list(fill="white"))+
  stat_compare_means(label.y = 11,method = "anova")+     # Global p-value
  stat_compare_means(comparisons = my_comparisons,label.y = c(7,9,8))

ggbarplot(PlantGrowth, x = "group", y = "weight", add = "mean_se",
          color = "group",fill="group", palette = "jco")+
  stat_compare_means(label.y = 11,method = "anova") +    # Global p-value
  stat_compare_means(comparisons = my_comparisons,label.y = c(7,9, 8),method="t.test")+
  scale_y_continuous(expand=c(0,0))

或者線圖

ggline(my_data, x = "group", y = "weight", 
       add = c("mean_se", "jitter"), 
       order = c("ctrl", "trt1", "trt2"),
       ylab = "Weight", xlab = "Treatment")+
      stat_compare_means(label.y = 11,method = "anova")+   # Global p-value
      stat_compare_means(comparisons = my_comparisons,label.y = c(7,9,8))

以上就是“R語言如何實現(xiàn)方差分析及其可視化”這篇文章的所有內(nèi)容，感謝各位的閱讀！相信大家閱讀完這篇文章都有很大的收獲，小編每天都會為大家更新不同的知識，如果還想學習更多的知識，請關(guān)注億速云行業(yè)資訊頻道。