R plotRresearch

R可视化重复测量的方差分析,自动两两比较

2023-05-06  本文已影响0人  欧阳松

参考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
上一篇下一篇

猜你喜欢

热点阅读