本篇推文来自于公众号读者的投稿，编辑排版由小明完成
1、单因素方差分析
1.1 加载R包
代码语言：javascript复制library(ggpubr)
library(rstatix)
library(tidyverse)
1.2 数据准备
这里用到的是R语言的内置数据集sample_n_by()函数很有用，能够分组随机抽样%>% 是管道符 是将前面的结果传输给后面的函数
data("PlantGrowth")
set.seed(1234)
PlantGrowth %>% sample_n_by(group, size = 1)
函数sample_n_by()加载和检查数据，按组显示随机的一行
显示分组变量的levels
代码语言：javascript复制levels(PlantGrowth$group)
单因素方差分析可以用来确定在三种条件下植物的平均生长是否显著不同。
1、3 统计
按组计算均值与标准差
代码语言：javascript复制PlantGrowth %>%
group_by(group) %>%
get_summary_stats(weight, type = "mean_sd")
1.4正态性假设
QQ图
代码语言：javascript复制model <- lm(weight ~ group, data = PlantGrowth)
ggqqplot(residuals(model))
image.png
Shapiro-Wilk test
代码语言：javascript复制shapiro_test(residuals(model))
在QQ图中，由于所有的点都近似地落在参考线上，我们可以假设正态分布。这个结论得到了Shapiro-Wilk test的支持。p值不显著(p=0.13>0.05)，因此我们可以假设正态性。
分组正态性检验
代码语言：javascript复制PlantGrowth %>%
group_by(group) %>%
shapiro_test(weight)
p > 0.05 假设成立
分组qq图
代码语言：javascript复制ggqqplot(PlantGrowth, "weight", facet.by = "group")
image.png
1.5方差同一性假设
残差拟合图可用于检验方差齐性。
代码语言：javascript复制plot(model, 1)
image.png
在上图中，残差与拟合值(每组的均值)之间没有明显的关系。我们可以假设方差齐性。
Levene’s test 检查方差齐性
代码语言：javascript复制PlantGrowth %>% levene_test(weight ~ group)
p>0.05, 没有显著性差异，假设通过。
1.6计算
代码语言：javascript复制res.aov <- PlantGrowth %>% anova_test(weight ~ group)
res.aov
p<0.05,表明各组之间有显著性差异
1.7多重比较
代码语言：javascript复制pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)
pwc
1.8带有p值的箱型图
代码语言：javascript复制pwc <- pwc %>% add_xy_position(x = "group")
ggboxplot(PlantGrowth, x = "group", y = "weight",color = "#488dcb")
stat_pvalue_manual(pwc, hide.ns = TRUE)
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)
image.png
2、双因素方差分析
代码语言：javascript复制library(ggpubr)
library(rstatix)
library(tidyverse)
2.1加载数据并按组检查任意一行
代码语言：javascript复制set.seed(123)
library(datarium)
data("jobsatisfaction", package = "datarium")
jobsatisfaction %>% sample_n_by(gender, education_level, size = 1)
2.2统计
均值与方差
代码语言：javascript复制jobsatisfaction %>%
group_by(gender, education_level) %>%
get_summary_stats(score, type = "mean_sd")
2.3可视化
代码语言：javascript复制bxp <- ggboxplot(
jobsatisfaction, x = "gender", y = "score",
color = "education_level", palette = "jco"
)
bxp
image.png
2.4正态性假设
建立线性模型
代码语言：javascript复制model1 <- lm(score ~ gender*education_level,
data = jobsatisfaction)
创建QQ图
代码语言：javascript复制ggqqplot(residuals(model1))
image.png
Shapiro-Wilk test
代码语言：javascript复制shapiro_test(residuals(model1))
假设通过
按组检验正态性
代码语言：javascript复制jobsatisfaction %>%
group_by(gender, education_level) %>%
shapiro_test(score)
p>0.05,假设通过
分组QQ图
代码语言：javascript复制ggqqplot(jobsatisfaction, "score", ggtheme = theme_bw())
facet_grid(gender ~ education_level)
image.png
2.5方差同一性假设
代码语言：javascript复制jobsatisfaction %>% levene_test(score ~ gender*education_level)
p > 0.05假设通过
2.6计算
代码语言：javascript复制res.aov1 <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov1
性别与教育程度对工作满意度得分的交互作用有统计学意义
2.7多重检验
简单主效应比较
代码语言：javascript复制model3 <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
group_by(gender) %>%
anova_test(score ~ education_level, error = model3)
受教育程度”对工作满意度的简单主效应在男性和女性中均有统计学意义。
两两比较
代码语言：javascript复制install.packages("emmeans")
library(emmeans)
pwc1 <- jobsatisfaction %>%
group_by(gender) %>%
emmeans_test(score ~ education_level, p.adjust.method = "bonferroni")
pwc1
各组男性和女性的工作满意度得分均有显著性差异(p < 0.05)
2.8可视化
代码语言：javascript复制pwc1 <- pwc1 %>% add_xy_position(x = "gender")
bxp
stat_pvalue_manual(pwc1)
labs(
subtitle = get_test_label(res.aov1, detailed = TRUE),
caption = get_pwc_label(pwc1)
)
image.png
原文链接：https://www.datanovia.com/en/checkout/order-received/
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
本文的完整代码可以在公众号后台回复 20210817 获得
