参考Mixed ANOVA in R: The Ultimate Guide - Datanovia
library(tidyverse)
library(ggpubr)
library(rstatix)
数据准备
set.seed(123)
data("anxiety", package = "datarium")
anxiety %>% sample_n_by(group, size = 1)
将短数据转换为长数据
# 将t1, t2 和 t3 列转换为长数据格式
# 将id和time列转换为因子形式
anxiety <- anxiety %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
# Inspect some random rows of the data by groups
set.seed(123)
anxiety %>% sample_n_by(group, time, size = 1)
数据统计
按照time和group分组，计算score的总结统计（均数和标准差）
anxiety %>%
group_by(time, group) %>%
get_summary_stats(score, type = "mean_sd")
## # A tibble: 9 x 6
## group time variable n mean sd
## <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 grp1 t1 score 15 17.1 1.63
## 2 grp2 t1 score 15 16.6 1.57
## 3 grp3 t1 score 15 17.0 1.32
## 4 grp1 t2 score 15 16.9 1.70
## 5 grp2 t2 score 15 16.5 1.70
## 6 grp3 t2 score 15 15.0 1.39
## # … with 3 more rows
可视化
bxp <- ggboxplot(
anxiety, x = "time", y = "score",
color = "group", palette = "jco"
)
bxp
boxplot
假设性检验
异常值
anxiety %>%
group_by(time, group) %>%
identify_outliers(score)
## [1] group time id score is.outlier is.extreme
## <0 rows> (or 0-length row.names)
正态性假设
anxiety %>%
group_by(time, group) %>%
shapiro_test(score)
## group time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 grp1 t1 score 0.964 0.769
## 2 grp2 t1 score 0.977 0.949
## 3 grp3 t1 score 0.954 0.588
## 4 grp1 t2 score 0.956 0.624
## 5 grp2 t2 score 0.935 0.328
## 6 grp3 t2 score 0.952 0.558
可视化
ggqqplot(anxiety, "score", ggtheme = theme_bw()) +
facet_grid(time ~ group)
qqplot
方差齐性假设
anxiety %>%
group_by(time) %>%
levene_test(score ~ group)
## # A tibble: 3 x 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 t1 2 42 0.176 0.839
## 2 t2 2 42 0.249 0.781
## 3 t3 2 42 0.335 0.717
计算
# Two-way mixed ANOVA test
res.aov <- anova_test(
data = anxiety, dv = score, wid = id,
between = group, within = time
)
get_anova_table(res.aov)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 group 2 42 4.35 1.90e-02 * 0.168
## 2 time 2 84 394.91 1.91e-43 * 0.179
## 3 group:time 4 84 110.19 1.38e-32 * 0.108
事后检验
# Effect of group at each time point
one.way <- anxiety %>%
group_by(time) %>%
anova_test(dv = score, wid = id, between = group) %>%
get_anova_table() %>%
adjust_pvalue(method = "bonferroni")
one.way
## # A tibble: 3 x 9
## time Effect DFn DFd F p `p<.05` ges p.adj
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 t1 group 2 42 0.365 0.696 "" 0.017 1
## 2 t2 group 2 42 5.84 0.006 * 0.218 0.018
## 3 t3 group 2 42 13.8 0.0000248 * 0.396 0.0000744
# Pairwise comparisons between group levels
pwc <- anxiety %>%
group_by(time) %>%
pairwise_t_test(score ~ group, p.adjust.method = "bonferroni")
pwc
## # A tibble: 9 x 10
## time .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 t1 score grp1 grp2 15 15 0.43 ns 1 ns
## 2 t1 score grp1 grp3 15 15 0.895 ns 1 ns
## 3 t1 score grp2 grp3 15 15 0.51 ns 1 ns
## 4 t2 score grp1 grp2 15 15 0.435 ns 1 ns
## 5 t2 score grp1 grp3 15 15 0.00212 ** 0.00636 **
## 6 t2 score grp2 grp3 15 15 0.0169 * 0.0507 ns
## # … with 3 more rows
可视化报告
# 可视化: boxplots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
pwc.filtered <- pwc %>% filter(time != "t1")
bxp +
stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)
image.png
