R海拾遗-单项重复测量方差分析

火星娃统计

发布于 2020-09-15 15:39:28

1.9K0

发布于 2020-09-15 15:39:28

文章被收录于专栏：火星娃统计火星娃统计火星娃统计

One-way repeated measures ANOVA

sunqi

2020/7/24

概述

重复测量资料总共包含下面三个部分，此次为第一个部分

One-way repeated measures ANOVA
two-way repeated measures ANOVA
three-way repeated measures ANOVA

假设

无明显的异常值使用identify_outliers()[rstatix package]
正态性：使用 Shapiro-Wilk normality test (shapiro_test() [rstatix]) 或者Shapiro-Wilk normality test (shapiro_test() [rstatix])
球形假设：通过anova_test()自动汇报 [rstatix package]

对于单项的方差分析，如果不满足上述假设，使用Friedman test进行，对于双向、三向的方差分析没有替代的非参数方法，只能通过装换数据

分析

需要的包

tidyverse：数据操作
ggpubr ：绘图
rstatix：管道符号
datarium：测试集数据

# 如果需要请安装
library(tidyverse)
library(ggpubr)
library(rstatix)

one-way

# 获得示例数据
data("selfesteem", package = "datarium")
head(selfesteem, 3)

## # A tibble: 3 x 4
##      id    t1    t2    t3
##   <int> <dbl> <dbl> <dbl>
## 1     1  4.01  5.18  7.11
## 2     2  2.56  6.91  6.31
## 3     3  3.24  4.44  9.78

#将列t1、t2和t3转换为长格式
#将id和time转换为因素变量

selfesteem <- selfesteem %>%
  gather(key = "time", value = "score", t1, t2, t3) %>%
  convert_as_factor(id, time)
head(selfesteem, 3)

## # A tibble: 3 x 3
##   id    time  score
##   <fct> <fct> <dbl>
## 1 1     t1     4.01
## 2 2     t1     2.56
## 3 3     t1     3.24

# 描述数据
selfesteem %>%
  group_by(time) %>%
  get_summary_stats(score, type = "mean_sd")

## # A tibble: 3 x 5
##   time  variable     n  mean    sd
##   <fct> <chr>    <dbl> <dbl> <dbl>
## 1 t1    score       10  3.14 0.552
## 2 t2    score       10  4.93 0.863
## 3 t3    score       10  7.64 1.14

# 可视化
bxp <- ggboxplot(selfesteem, x = "time", y = "score", add = "point")
bxp

# 检验假设
# 异常值检验
selfesteem %>%
  group_by(time) %>%
  identify_outliers(score)

## # A tibble: 2 x 5
##   time  id    score is.outlier is.extreme
##   <fct> <fct> <dbl> <lgl>      <lgl>
## 1 t1    6      2.05 TRUE       FALSE
## 2 t2    2      6.91 TRUE       FALSE

# 没有极端的异常值

# 正态假设
selfesteem %>%
  group_by(time) %>%
  shapiro_test(score)

## # A tibble: 3 x 4
##   time  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 t1    score        0.967 0.859
## 2 t2    score        0.876 0.117
## 3 t3    score        0.923 0.380

# 注意，如果你的样本量大于50，建议使用QQ图
# 因为在较大的样本量下，Shapiro-Wilk测试变得非常敏感，即使是一个小的偏离正常值

# qq图
ggqqplot(selfesteem, "score", facet.by = "time")

# 球形假设
# 对于违反球形假设的数据，Greenhouse-Geisser sphericity 校正自动进行
res.aov <- anova_test(data = selfesteem, dv = score, wid = id, within = time)
get_anova_table(res.aov)

## ANOVA Table (type III tests)
##
##   Effect DFn DFd      F        p p<.05   ges
## 1   time   2  18 55.469 2.01e-08     * 0.829

ges：是广义效应量，受试体内因素影响的变异量

# 事后检验
pwc <- selfesteem %>%
  pairwise_t_test(
    score ~ time, paired = TRUE,
    p.adjust.method = "bonferroni"
    )
pwc

## # A tibble: 3 x 10
##   .y.   group1 group2    n1    n2 statistic    df         p   p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>   <dbl> <chr>
## 1 score t1     t2        10    10     -4.97     9   7.72e-4 2.00e-3 **
## 2 score t1     t3        10    10    -13.2      9   3.34e-7 1.00e-6 ****
## 3 score t2     t3        10    10     -4.87     9   8.86e-4 3.00e-3 **

# 可视化

# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
bxp +
  stat_pvalue_manual(pwc) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )