R语言单、双因素方差分析及结果可视化的简单小例子

用户7010445

发布于 2021-08-31 17:47:50

6.2K10

代码可运行

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

运行总次数：0

代码可运行

本篇推文来自于公众号读者的投稿，编辑排版由小明完成

1、单因素方差分析

1.1 加载R包

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

levels(PlantGrowth$group)

单因素方差分析可以用来确定在三种条件下植物的平均生长是否显著不同。

1、3 统计

按组计算均值与标准差

PlantGrowth %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "mean_sd")

1.4正态性假设

QQ图

model  <- lm(weight ~ group, data = PlantGrowth)
ggqqplot(residuals(model))

image.png

Shapiro-Wilk test

shapiro_test(residuals(model))

在QQ图中，由于所有的点都近似地落在参考线上，我们可以假设正态分布。这个结论得到了Shapiro-Wilk test的支持。p值不显著(p=0.13>0.05)，因此我们可以假设正态性。

分组正态性检验

PlantGrowth %>%
  group_by(group) %>%
  shapiro_test(weight)

p > 0.05 假设成立

分组qq图

ggqqplot(PlantGrowth, "weight", facet.by = "group")

image.png

1.5方差同一性假设

残差拟合图可用于检验方差齐性。

plot(model, 1)

image.png

在上图中，残差与拟合值(每组的均值)之间没有明显的关系。我们可以假设方差齐性。

Levene’s test 检查方差齐性

PlantGrowth %>% levene_test(weight ~ group)

p>0.05, 没有显著性差异，假设通过。

1.6计算

res.aov <- PlantGrowth %>% anova_test(weight ~ group)
res.aov

p<0.05,表明各组之间有显著性差异

1.7多重比较

pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)
pwc

1.8带有p值的箱型图

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、双因素方差分析

library(ggpubr)
library(rstatix)
library(tidyverse)

2.1加载数据并按组检查任意一行

set.seed(123)
library(datarium)
data("jobsatisfaction", package = "datarium")
jobsatisfaction %>% sample_n_by(gender, education_level, size = 1)

2.2统计

均值与方差

jobsatisfaction %>%
  group_by(gender, education_level) %>%
  get_summary_stats(score, type = "mean_sd")

2.3可视化

bxp <- ggboxplot(
  jobsatisfaction, x = "gender", y = "score",
  color = "education_level", palette = "jco"
)
bxp

image.png

2.4正态性假设

建立线性模型

model1  <- lm(score ~ gender*education_level,
             data = jobsatisfaction)

创建QQ图

ggqqplot(residuals(model1))

image.png

Shapiro-Wilk test

shapiro_test(residuals(model1))

假设通过

按组检验正态性

jobsatisfaction %>%
  group_by(gender, education_level) %>%
  shapiro_test(score)

p>0.05,假设通过

分组QQ图

ggqqplot(jobsatisfaction, "score", ggtheme = theme_bw()) +
  facet_grid(gender ~ education_level)

image.png

2.5方差同一性假设

jobsatisfaction %>% levene_test(score ~ gender*education_level)

p > 0.05假设通过

2.6计算

res.aov1 <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov1

性别与教育程度对工作满意度得分的交互作用有统计学意义

2.7多重检验

简单主效应比较

model3 <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model3)

受教育程度”对工作满意度的简单主效应在男性和女性中均有统计学意义。

两两比较

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可视化

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 获得

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-08-17，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

r 语言

本文分享自小明的数据分析笔记本微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据分析

r 语言

登录后参与评论

0 条评论

热度