r library car_R绘图应用实例：单因素方差分析ANOVA及绘图-CSDN博客

本文主要是利用日常实验数据，尝试用R进行单因素方差分析并绘制柱形图。

ANOVA原理参考：单因素方差分析(One-way Anova)实验数据：
在随机划分的试验田中，施加三种复合肥(B,C,D)，饲料填充物做空白对照(A)，一段时间后，测定试验田内植株高度，比较数据有无不同，差异性是否具有统计学意义。
计算结果录入Excel，如下：

图片.png

(一)：数据准备

rm(list = ls())
# prepare the data
data                    B=c(68,55,49,62,70,59,63),
                   C=c(78,76,72,81,76,71,83),
                   D=c(85,65,81,98,75,92,79))
head(data)

# convert data from a wide format to a long format.
library(reshape2)
data_long head(data_long)

(二)：方差分析(条件适用检验及方差分析)

# 单因素方差分析(One-way Anova)前提假设的检验：
# 1. 正态性检验(对4个水平下的每组数据都做一次正态检验)
shapiro.test(data$A) #或shapiro.test(data_long$value[1:7])
shapiro.test(data$B) #或shapiro.test(data_long$value[8:14])
shapiro.test(data$C) #或shapiro.test(data_long$value[15:21])
shapiro.test(data$D) #或shapiro.test(data_long$value[22:28])
############################################################################
#split:拆分为list数据集;                                                    #
#lapply:对数据集进行shapiro.test检验循环;                                   #
#unlist():向量的原样返回                                                    #
group_data group_data                                                                 #
unlist(lapply(group_data, function(x){                                     #
  shapiro.test(x)$p.value                                                  #
}))                                                                        #
############################################################################
# 使用Q-Q图来检验正态性
library(car)
qqPlot(group_data[[1]]) #或qqPlot(data$A)
qqPlot(group_data[[2]]) #或qqPlot(data$A)
qqPlot(group_data[[3]]) #或qqPlot(data$A)
qqPlot(group_data[[4]]) #或qqPlot(data$D)

# 2. 方差齐性检验，使用car包的leveneTest()
leveneTest(value~variable, data = data_long)

# 3. 离群点检验：
outlierTest(lm(value~variable, data = data_long))

Q-Q图示例

# 4. 单因素方差分析ANOVA
aov1 summary(aov1)
            Df Sum Sq Mean Sq F value   Pr(>F)    
variable     3   4049  1349.7    18.2 2.25e-06 ***
Residuals   24   1780    74.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Show the means
model.tables(aov1, "means")

# Tukey HSD post-hoc test
TukeyHSD(aov1)
  Tukey multiple comparisons of means
    95% family-wise confidence level
Fit: aov(formula = value ~ variable, data = data_long)
$variable
         diff       lwr      upr     p adj
B-A  8.714286 -3.984450 21.41302 0.2574525
C-A 24.571429 11.872693 37.27016 0.0000986
D-A 30.000000 17.301265 42.69874 0.0000055
C-B 15.857143  3.158407 28.55588 0.0106071
D-B 21.285714  8.586979 33.98445 0.0005874
D-C  5.428571 -7.270164 18.12731 0.6453377

# Alternative: Multiple comparisons using multcomp package
library(multcomp)
e summary(e)
     Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrast
Fit: aov(formula = value ~ variable, data = data_long)
Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
B - A == 0    8.714      4.603   1.893   0.2573    
C - A == 0   24.571      4.603   5.338   <0.001 ***
D - A == 0   30.000      4.603   6.517   <0.001 ***
C - B == 0   15.857      4.603   3.445   0.0106 *  
D - B == 0   21.286      4.603   4.624   <0.001 ***
D - C == 0    5.429      4.603   1.179   0.6453    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

(三)：统计分析

# summarySE 函数提供了标准差、标准误以及95%的置信区间
library(Rmisc)
data_long_count                                groupvars= "variable")
data_long_count

(四)：绘图

绘图应用1

# 绘制带有显著性标记的条形图
marker library(ggplot2)
ggplot(data_long_count, aes(x=variable, y=value,fill=variable)) + 
  geom_bar(stat="identity", color="black",size=.3) +  
  geom_errorbar(aes(ymin=value-se, ymax=value+se),
                size=.3, width=.2) +
  geom_text(aes(y = value +  1.5 * se, label = marker), 
            size = 5, fontface = "bold") +
  scale_color_brewer("Set1") +
  ggtitle("Effects of Three Compound Fertilizers on Plant Height") +
  xlab("Different compound fertilizer treatments") +
  ylab("Plant height")

Rplot04.png

绘图应用2

# 最终绘图效果
library(ggsignif) 
library(ggplot2)
ggplot(data_long_count, aes(x=variable, y=value,fill=variable)) + 
  geom_bar(stat="identity", color="black",size=.3) + 
  # 添加误差棒
  geom_errorbar(aes(ymin=value-se, ymax=value+se), size=.3, width=.2) +
  # 手动添加显著性
  geom_signif(annotations = c("***","***"), y_position = c(83, 90),
              xmin = c(1, 1), xmax = c(3, 4), 
              tip_length = c(c(0.65, 0.05),c(0.7, 0.05)), vjust = 0 ) +
  # 设置色板
  scale_color_brewer(palette = "Set1") +
  # 更改主标题
  ggtitle("Effects of Three Compound Fertilizers on Plant Height") +
  xlab("Different compound fertilizer treatments") +
  ylab("Plant height")  +
  # 更改轴标签的外观
  theme(plot.title = element_text(color="red", size=14, face="bold.italic"),
        axis.title.x = element_text(color="blue", size=14, face="bold"),
        axis.title.y = element_text(color="#993333", size=14, face="bold")) +
  # 更改图例标签
  labs(fill="Treated")