8  重复测量方差分析

Modified

November 20, 2024

重复测量方差分析用于分析重复测量数据中组间和组内自变量对因变量的影响。模型形式为：

\[ Y_{ijk}=\mu +\alpha_i +\beta_j +(\alpha\beta)_{ij}+S_k+(\beta S)_{jk}+\epsilon_{ijk} \] 其中，\(Y_{ijk}\) 是第i 个组、第j 个时间点或条件下、第k 个受试者的测量值，\(\alpha_i\) 是组间效应，\(\beta_j\) 时间效应或条件效应，\((\alpha\beta)_{ij}\) 是组效应与时间效应的交互效应，\(S_k\) 是第 k 个受试者的效应（通常被视为随机效应），\((\beta S)_{jk}\) 是时间效应与受试者效应的交互作用（也通常被视为随机效应）， \(\epsilon_{ijk}\) 是误差项。

8.1 假设

  • 独立性，即不同受试者之间的观测值是独立的。

  • 正态性，观测值在每个时间点或条件下的正态分布。

  • 同一受试者的重复测量之间的相关性一致性。

  • 各组内的方差齐性（sphericity），即组内相关性结构是已知的，通常假设为复合对称（compound symmetry）。协方差矩阵的球形检验（W=1）

Note

  • 复合对称假设：RM-ANOVA假设误差项具有复合对称性，即组内相关性和方差齐性。如果复合对称假设不成立，结果可能不准确，可以使用更为灵活的线性混合模型（LMM）或广义估计方程（GEE）进行分析。

  • 检验方法：RM-ANOVA的主要检验方法包括F检验，用于检验组效应、时间效应和交互作用是否显著。

8.2 单因素重复测量方差分析

8.2.1 公式

\[ SS_T=\sum_{j=1}^{t}\sum_{i=1}^{n}X_{ij}^2-\frac{(\sum_{j=1}^{t}\sum_{i=1}^{n}X_{ij})^2}{nt}=SS_{受试者间}+SS_{不同时间点}+SS_{E\ \ 随机误差}，\nu_T=nt-1 \]

\[ SS_{受试者间}=\frac{\sum_{i}(\sum_{j}X_{ij})^2}{t}-C，\nu_{受试者}=n-1 \]反映了在受试者间的变异。

\[ SS_{不同时间点}=\frac{\sum_{j}(\sum_{i}X_{ij})^2}{n}-C ，\nu_{时间点}=t-1 \]反映了在每个受试者内不同时间点的重复测量变异。

\[ SS_E=SS_T-SS_{受试者间}-SS_{不同时间点}，\nu_E=(n-1)(t-1) \]

\(H_0:\mu_1=\mu_2=...=\mu_t\)，检验统计量

\[ F=\frac{MS_{不同时间点}}{MS_E}=\frac{SS_{不同时间点}/(t-1)}{SS_E/((n-1)(t-1))} \]

Show the code
ALT <- tribble(
    ~id,~zero,~ten,~twenty,~thirty,
    1,186,122,134,110,
    2,345,312,268,176,
    3,98,84,52,61,
    4,288,98,91,85,
    5,176,86,130,99,
    6,210,188,143,120,
    7,271,322,86,65,
    8,415,332,265,186,
    9,171,126,130,135,
    10,243,330,95,64,
)
ALT_long <-
    ALT |> pivot_longer(cols = -1,
                       names_to = "week",
                       values_to = "ALT") |> 
    mutate(week = factor(week, levels = c("zero","ten", "twenty", "thirty")))


# 假设检验
ALT_long %>% 
    group_by(week) %>% 
    rstatix::shapiro_test(ALT)
#> # A tibble: 4 × 4
#>   week   variable statistic      p
#>   <fct>  <chr>        <dbl>  <dbl>
#> 1 zero   ALT          0.972 0.906 
#> 2 ten    ALT          0.796 0.0131
#> 3 twenty ALT          0.834 0.0378
#> 4 thirty ALT          0.907 0.262
Show the code
# 总变异
xij_sum <- sum(ALT_long$ALT)
xij_square_sum <- sum(ALT_long$ALT^2)
C <- xij_sum^2/40
SS_T <- xij_square_sum-C

# 受试者间  
xi._sum <- rowSums(ALT[,-1])
xi._sum
#>  [1]  552 1101  295  562  491  661  744 1198  562  732
SS_B <- sum(xi._sum^2)/4-C #行和平方和
MS_B <- SS_B/9

# 不同时间点
x.j_sum <- colSums(ALT[,-1])
x.j_sum
#>   zero    ten twenty thirty 
#>   2403   2000   1394   1101
SS_W <- sum(x.j_sum^2)/10-C #列和平方和
MS_W <- SS_W/3



SS_E <- SS_T-SS_B-SS_W
MS_E <- SS_E/(9*3)

F_stat <- MS_W/MS_E
F_stat
#> [1] 11.13855

v1 = 4-1
v2 = (4-1)*(10-1)

pf(F_stat, df1 = v1, df2 = v2, lower.tail = F)
#> [1] 6.170758e-05

cat("F统计量：", F_stat, "\n",
    "p值：",pf(F_stat, df1 = v1, df2 = v2, lower.tail = F))
#> F统计量： 11.13855 
#>  p值： 6.170758e-05

8.2.2 nlme::lme()

Show the code
library(nlme)
ALT_lme <- lme(fixed=ALT ~ week, random = ~1| id/week, data = ALT_long)

ALT_aov <- anova.lme(ALT_lme)
ALT_aov 
#>             numDF denDF  F-value p-value
#> (Intercept)     1    27 62.98194  <.0001
#> week            3    27 11.13855   1e-04
ALT_aov$`p-value`
#> [1] 1.569667e-08 6.170758e-05

事后检验

Show the code
#成对比较
library(emmeans)
emmeans(ALT_lme, pairwise ~week, adjust = "bonferroni") 
#> $emmeans
#>  week   emmean   SE df lower.CL upper.CL
#>  zero      240 26.5  9    180.3      300
#>  ten       200 26.5  9    140.0      260
#>  twenty    139 26.5  9     79.4      199
#>  thirty    110 26.5  9     50.1      170
#> 
#> Degrees-of-freedom method: containment 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#>  contrast        estimate   SE df t.ratio p.value
#>  zero - ten          40.3 24.9 27   1.620  0.7014
#>  zero - twenty      100.9 24.9 27   4.055  0.0023
#>  zero - thirty      130.2 24.9 27   5.233  0.0001
#>  ten - twenty        60.6 24.9 27   2.436  0.1305
#>  ten - thirty        89.9 24.9 27   3.613  0.0073
#>  twenty - thirty     29.3 24.9 27   1.178  1.0000
#> 
#> Degrees-of-freedom method: containment 
#> P value adjustment: bonferroni method for 6 tests
emmeans(ALT_lme,  ~week) %>% pairs(., adjust = "bonferroni")
#>  contrast        estimate   SE df t.ratio p.value
#>  zero - ten          40.3 24.9 27   1.620  0.7014
#>  zero - twenty      100.9 24.9 27   4.055  0.0023
#>  zero - thirty      130.2 24.9 27   5.233  0.0001
#>  ten - twenty        60.6 24.9 27   2.436  0.1305
#>  ten - thirty        89.9 24.9 27   3.613  0.0073
#>  twenty - thirty     29.3 24.9 27   1.178  1.0000
#> 
#> Degrees-of-freedom method: containment 
#> P value adjustment: bonferroni method for 6 tests
Show the code
library(multcomp)
ALT_posthoc <-glht(ALT_lme, linfct = mcp(week =c("zero - ten = 0 ",
                                               "twenty - zero == 0",
                                               "thirty - zero = 0 ")))
summary(ALT_posthoc) %>% broom::tidy()
#> # A tibble: 3 × 7
#>   term  contrast      null.value estimate std.error statistic adj.p.value
#>   <chr> <chr>              <dbl>    <dbl>     <dbl>     <dbl>       <dbl>
#> 1 week  zero - ten             0     40.3      24.9      1.62 0.248      
#> 2 week  twenty - zero          0   -101.       24.9     -4.06 0.000167   
#> 3 week  thirty - zero          0   -130.       24.9     -5.23 0.000000418

8.3 含组间因子的重复测量方差分析

8.3.1 公式

\[ \begin{aligned} SS_总&= SS_{受试对象间}+SS_{受试对象内} \\ &=（SS_{处理方法}+SS_{个体间差异}）+（SS_{时间}+SS_{处理与时间交互}+SS_{个体内差异}） \end{aligned} \]

https://personality-project.org/r/r.guide/r.anova.html#oneway

8.3.2 案例

孙振球（医学统计学第四版 人卫版，例12-3，P193）

Show the code
# 医学统计学 人卫版 第7版
SBP <- read_tsv("data/麻醉诱导时相.tsv")
SBP
#> # A tibble: 15 × 7
#>       id group    t0    t1    t2    t3    t4
#>    <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1     1 A       120   108   112   120   117
#>  2     2 A       118   109   115   126   123
#>  3     3 A       119   112   119   124   118
#>  4     4 A       121   112   119   126   120
#>  5     5 A       127   121   127   133   126
#>  6     6 B       121   120   118   131   137
#>  7     7 B       122   121   119   129   133
#>  8     8 B       128   129   126   135   142
#>  9     9 B       117   115   111   123   131
#> 10    10 B       118   114   116   123   133
#> 11    11 C       131   119   118   135   129
#> 12    12 C       129   128   121   148   132
#> 13    13 C       123   123   120   143   136
#> 14    14 C       123   121   116   145   126
#> 15    15 C       125   124   118   142   130
Show the code
SBP_long <- SBP |> pivot_longer(cols = starts_with("t"),
                              names_to = "time",
                              values_to = "SBP") |> 
    mutate(
        id=factor(id),
        group=factor(group),
        time=factor(time)
    )

SBP_long %>% head()
#> # A tibble: 6 × 4
#>   id    group time    SBP
#>   <fct> <fct> <fct> <dbl>
#> 1 1     A     t0      120
#> 2 1     A     t1      108
#> 3 1     A     t2      112
#> 4 1     A     t3      120
#> 5 1     A     t4      117
#> 6 2     A     t0      118
Show the code
# 正态性检验
library(rstatix)
SBP_long %>% group_by(time,group) %>% 
    shapiro_test(SBP)
#> # A tibble: 15 × 5
#>    group time  variable statistic     p
#>    <fct> <fct> <chr>        <dbl> <dbl>
#>  1 A     t0    SBP          0.836 0.154
#>  2 B     t0    SBP          0.916 0.503
#>  3 C     t0    SBP          0.867 0.254
#>  4 A     t1    SBP          0.834 0.148
#>  5 B     t1    SBP          0.913 0.485
#>  6 C     t1    SBP          0.978 0.921
#>  7 A     t2    SBP          0.940 0.666
#>  8 B     t2    SBP          0.969 0.869
#>  9 C     t2    SBP          0.953 0.758
#> 10 A     t3    SBP          0.937 0.647
#> 11 B     t3    SBP          0.902 0.421
#> 12 C     t3    SBP          0.941 0.672
#> 13 A     t4    SBP          0.943 0.687
#> 14 B     t4    SBP          0.892 0.367
#> 15 C     t4    SBP          0.984 0.955

8.3.2.1 nlme::lme()

Show the code
# 第一种方法
library(nlme)
SBP_lme <- lme(SBP ~ group*time, random = ~ 1 |id/time, data = SBP_long)
# summary(SBP_lme)
anova(SBP_lme)
#>             numDF denDF   F-value p-value
#> (Intercept)     1    48 14649.223  <.0001
#> group           2    12     5.783  0.0174
#> time            4    48   106.558  <.0001
#> group:time      8    48    19.101  <.0001
# 看不到球形检验


# 事后检验 成对比较
library(emmeans)
# 如果交互效应不显著

## 时间点事后检验
SBP_long %>%
  pairwise_t_test(SBP ~ time, p.adjust.method = "none", detailed = TRUE)
#> # A tibble: 10 × 10
#>    .y.   group1 group2    n1    n2        p method   p.adj p.signif p.adj.signif
#>  * <chr> <chr>  <chr>  <int> <int>    <dbl> <chr>    <dbl> <chr>    <chr>       
#>  1 SBP   t0     t1        15    15  6.81e-2 T-test 6.81e-2 ns       ns          
#>  2 SBP   t0     t2        15    15  6.41e-2 T-test 6.41e-2 ns       ns          
#>  3 SBP   t1     t2        15    15  9.78e-1 T-test 9.78e-1 ns       ns          
#>  4 SBP   t0     t3        15    15  1.78e-4 T-test 1.78e-4 ***      ***         
#>  5 SBP   t1     t3        15    15  1.67e-7 T-test 1.67e-7 ****     ****        
#>  6 SBP   t2     t3        15    15  1.49e-7 T-test 1.49e-7 ****     ****        
#>  7 SBP   t0     t4        15    15  1.28e-2 T-test 1.28e-2 *        *           
#>  8 SBP   t1     t4        15    15  3.68e-5 T-test 3.68e-5 ****     ****        
#>  9 SBP   t2     t4        15    15  3.32e-5 T-test 3.32e-5 ****     ****        
#> 10 SBP   t3     t4        15    15  1.65e-1 T-test 1.65e-1 ns       ns

emmeans(SBP_lme,pairwise ~time, adjust = "none")  
#> $emmeans
#>  time emmean   SE df lower.CL upper.CL
#>  t0      123 1.16 12      120      125
#>  t1      118 1.16 12      116      121
#>  t2      118 1.16 12      116      121
#>  t3      132 1.16 12      130      135
#>  t4      129 1.16 12      126      131
#> 
#> Results are averaged over the levels of: group 
#> Degrees-of-freedom method: containment 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#>  contrast estimate    SE df t.ratio p.value
#>  t0 - t1    4.4000 0.855 48   5.147  <.0001
#>  t0 - t2    4.4667 0.855 48   5.225  <.0001
#>  t0 - t3   -9.4000 0.855 48 -10.995  <.0001
#>  t0 - t4   -6.0667 0.855 48  -7.096  <.0001
#>  t1 - t2    0.0667 0.855 48   0.078  0.9382
#>  t1 - t3  -13.8000 0.855 48 -16.142  <.0001
#>  t1 - t4  -10.4667 0.855 48 -12.243  <.0001
#>  t2 - t3  -13.8667 0.855 48 -16.220  <.0001
#>  t2 - t4  -10.5333 0.855 48 -12.321  <.0001
#>  t3 - t4    3.3333 0.855 48   3.899  0.0003
#> 
#> Results are averaged over the levels of: group 
#> Degrees-of-freedom method: containment

## 组间事后检验
emmeans(SBP_lme,pairwise ~group, adjust = "none")  # 与spss一样
#> $emmeans
#>  group emmean   SE df lower.CL upper.CL
#>  A        120 1.78 14      116      123
#>  B        124 1.78 12      121      128
#>  C        128 1.78 12      124      132
#> 
#> Results are averaged over the levels of: time 
#> Degrees-of-freedom method: containment 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#>  contrast estimate   SE df t.ratio p.value
#>  A - B       -4.80 2.51 12  -1.911  0.0802
#>  A - C       -8.52 2.51 12  -3.392  0.0054
#>  B - C       -3.72 2.51 12  -1.481  0.1644
#> 
#> Results are averaged over the levels of: time 
#> Degrees-of-freedom method: containment

# 此处交互效应显著
# 1. 简单组别效应 与SPSS一样
SBP_long %>% group_by(time) %>% anova_test(SBP~group)
#> # A tibble: 5 × 8
#>   time  Effect   DFn   DFd      F        p `p<.05`   ges
#> * <fct> <chr>  <dbl> <dbl>  <dbl>    <dbl> <chr>   <dbl>
#> 1 t0    group      2    12  2.93  0.092    ""      0.328
#> 2 t1    group      2    12  6.03  0.015    "*"     0.501
#> 3 t2    group      2    12  0.022 0.979    ""      0.004
#> 4 t3    group      2    12 17.0   0.000312 "*"     0.74 
#> 5 t4    group      2    12 17.4   0.000286 "*"     0.743

# 2. 各时间点处理组成对比较  用的是未调整的p值，无法观察到标准误
SBP_long %>%
  group_by(time) %>%
  pairwise_t_test(SBP ~ group, 
                  p.adjust.method = "bonferroni", detailed = TRUE)
#> # A tibble: 15 × 11
#>    time  .y.   group1 group2    n1    n2         p method    p.adj p.signif
#>  * <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <chr>     <dbl> <chr>   
#>  1 t0    SBP   A      B          5     5 0.936     T-test 1        ns      
#>  2 t0    SBP   A      C          5     5 0.0539    T-test 0.162    ns      
#>  3 t0    SBP   B      C          5     5 0.0623    T-test 0.187    ns      
#>  4 t1    SBP   A      B          5     5 0.0358    T-test 0.107    *       
#>  5 t1    SBP   A      C          5     5 0.00541   T-test 0.0162   **      
#>  6 t1    SBP   B      C          5     5 0.327     T-test 0.981    ns      
#>  7 t2    SBP   A      B          5     5 0.894     T-test 1        ns      
#>  8 t2    SBP   A      C          5     5 0.947     T-test 1        ns      
#>  9 t2    SBP   B      C          5     5 0.842     T-test 1        ns      
#> 10 t3    SBP   A      B          5     5 0.456     T-test 1        ns      
#> 11 t3    SBP   A      C          5     5 0.000161  T-test 0.000483 ***     
#> 12 t3    SBP   B      C          5     5 0.000585  T-test 0.00175  ***     
#> 13 t4    SBP   A      B          5     5 0.0000887 T-test 0.000266 ****    
#> 14 t4    SBP   A      C          5     5 0.00201   T-test 0.00602  **      
#> 15 t4    SBP   B      C          5     5 0.0901    T-test 0.27     ns      
#> # ℹ 1 more variable: p.adj.signif <chr>

8.3.2.2 rstatix::anova_test() 看球形检验 W统计量

Show the code
library(rstatix)
aov_SBP <- anova_test(data = SBP_long,
           dv = SBP,
           wid = id,
           within = time,
           between = group,
           type = "3"
           )

aov_SBP
#> ANOVA Table (type III tests)
#> 
#> $ANOVA
#>       Effect DFn DFd       F        p p<.05   ges
#> 1      group   2  12   5.783 1.70e-02     * 0.430
#> 2       time   4  48 106.558 3.02e-23     * 0.659
#> 3 group:time   8  48  19.101 1.62e-12     * 0.409
#> 
#> $`Mauchly's Test for Sphericity`
#>       Effect     W     p p<.05
#> 1       time 0.293 0.178      
#> 2 group:time 0.293 0.178      
#> 
#> $`Sphericity Corrections`
#>       Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
#> 1       time 0.679 2.71, 32.58 1.87e-16         * 0.896 3.59, 43.03 4.65e-21
#> 2 group:time 0.679 5.43, 32.58 4.26e-09         * 0.896 7.17, 43.03 2.04e-11
#>   p[HF]<.05
#> 1         *
#> 2         *

# 进行事后检验，若交互效应不显著

SBP_long %>%
  group_by(time) %>%
  pairwise_t_test(SBP ~ group, p.adjust.method = "bonferroni", detailed = TRUE)
#> # A tibble: 15 × 11
#>    time  .y.   group1 group2    n1    n2         p method    p.adj p.signif
#>  * <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <chr>     <dbl> <chr>   
#>  1 t0    SBP   A      B          5     5 0.936     T-test 1        ns      
#>  2 t0    SBP   A      C          5     5 0.0539    T-test 0.162    ns      
#>  3 t0    SBP   B      C          5     5 0.0623    T-test 0.187    ns      
#>  4 t1    SBP   A      B          5     5 0.0358    T-test 0.107    *       
#>  5 t1    SBP   A      C          5     5 0.00541   T-test 0.0162   **      
#>  6 t1    SBP   B      C          5     5 0.327     T-test 0.981    ns      
#>  7 t2    SBP   A      B          5     5 0.894     T-test 1        ns      
#>  8 t2    SBP   A      C          5     5 0.947     T-test 1        ns      
#>  9 t2    SBP   B      C          5     5 0.842     T-test 1        ns      
#> 10 t3    SBP   A      B          5     5 0.456     T-test 1        ns      
#> 11 t3    SBP   A      C          5     5 0.000161  T-test 0.000483 ***     
#> 12 t3    SBP   B      C          5     5 0.000585  T-test 0.00175  ***     
#> 13 t4    SBP   A      B          5     5 0.0000887 T-test 0.000266 ****    
#> 14 t4    SBP   A      C          5     5 0.00201   T-test 0.00602  **      
#> 15 t4    SBP   B      C          5     5 0.0901    T-test 0.27     ns      
#> # ℹ 1 more variable: p.adj.signif <chr>

8.3.2.3 afex::aov_*()看各时间点处理组成对比较和各处理组时间点成对比较

https://cran.r-project.org/web/packages/afex/vignettes/afex_anova_example.html

Show the code
library("afex")     
library("emmeans") 
# aov_car(SBP ~ group + Error(id/time), data = SBP_long)
# aov_4(SBP ~ group + (time|id), data = SBP_long)

SBP_aovez <- aov_ez(id = "id",
                    dv =  "SBP",
                    data =  SBP_long, 
                    between ="group" , 
                    within = c("time")
                    )
SBP_aovez
#> Anova Table (Type 3 tests)
#> 
#> Response: SBP
#>       Effect          df   MSE          F  ges p.value
#> 1      group       2, 12 78.87     5.78 * .430    .017
#> 2       time 2.71, 32.58  8.08 106.56 *** .659   <.001
#> 3 group:time 5.43, 32.58  8.08  19.10 *** .409   <.001
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
#> 
#> Sphericity correction method: GG

SBP_aovez %>% summary()
#> 
#> Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
#> 
#>              Sum Sq num Df Error SS den Df    F value    Pr(>F)    
#> (Intercept) 1155433      1   946.48     12 14649.2234 < 2.2e-16 ***
#> group           912      2   946.48     12     5.7829   0.01743 *  
#> time           2336      4   263.12     48   106.5576 < 2.2e-16 ***
#> group:time      838      8   263.12     48    19.1006 1.621e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> 
#> Mauchly Tests for Sphericity
#> 
#>            Test statistic p-value
#> time              0.29307 0.17766
#> group:time        0.29307 0.17766
#> 
#> 
#> Greenhouse-Geisser and Huynh-Feldt Corrections
#>  for Departure from Sphericity
#> 
#>             GG eps Pr(>F[GG])    
#> time       0.67869  < 2.2e-16 ***
#> group:time 0.67869  4.263e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>              HF eps   Pr(>F[HF])
#> time       0.896371 4.649080e-21
#> group:time 0.896371 2.041727e-11

# 多变量检验
summary(SBP_aovez$Anova)
#> 
#> Type III Repeated Measures MANOVA Tests:
#> 
#> ------------------------------------------
#>  
#> Term: (Intercept) 
#> 
#>  Response transformation matrix:
#>    (Intercept)
#> t0           1
#> t1           1
#> t2           1
#> t3           1
#> t4           1
#> 
#> Sum of squares and products for the hypothesis:
#>             (Intercept)
#> (Intercept)     5777165
#> 
#> Multivariate Tests: (Intercept)
#>                  Df test stat approx F num Df den Df     Pr(>F)    
#> Pillai            1    0.9992 14649.22      1     12 < 2.22e-16 ***
#> Wilks             1    0.0008 14649.22      1     12 < 2.22e-16 ***
#> Hotelling-Lawley  1 1220.7686 14649.22      1     12 < 2.22e-16 ***
#> Roy               1 1220.7686 14649.22      1     12 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------
#>  
#> Term: group 
#> 
#>  Response transformation matrix:
#>    (Intercept)
#> t0           1
#> t1           1
#> t2           1
#> t3           1
#> t4           1
#> 
#> Sum of squares and products for the hypothesis:
#>             (Intercept)
#> (Intercept)      4561.2
#> 
#> Multivariate Tests: group
#>                  Df test stat approx F num Df den Df   Pr(>F)  
#> Pillai            2 0.4907894 5.782943      2     12 0.017434 *
#> Wilks             2 0.5092106 5.782943      2     12 0.017434 *
#> Hotelling-Lawley  2 0.9638239 5.782943      2     12 0.017434 *
#> Roy               2 0.9638239 5.782943      2     12 0.017434 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------
#>  
#> Term: time 
#> 
#>  Response transformation matrix:
#>    time1 time2 time3 time4
#> t0     1     0     0     0
#> t1     0     1     0     0
#> t2     0     0     1     0
#> t3     0     0     0     1
#> t4    -1    -1    -1    -1
#> 
#> Sum of squares and products for the hypothesis:
#>           time1     time2     time3     time4
#> time1  552.0667  952.4667  958.5333 -303.3333
#> time2  952.4667 1643.2667 1653.7333 -523.3333
#> time3  958.5333 1653.7333 1664.2667 -526.6667
#> time4 -303.3333 -523.3333 -526.6667  166.6667
#> 
#> Multivariate Tests: time
#>                  Df test stat approx F num Df den Df     Pr(>F)    
#> Pillai            1   0.98255 126.6592      4      9 6.6475e-08 ***
#> Wilks             1   0.01745 126.6592      4      9 6.6475e-08 ***
#> Hotelling-Lawley  1  56.29299 126.6592      4      9 6.6475e-08 ***
#> Roy               1  56.29299 126.6592      4      9 6.6475e-08 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------
#>  
#> Term: group:time 
#> 
#>  Response transformation matrix:
#>    time1 time2 time3 time4
#> t0     1     0     0     0
#> t1     0     1     0     0
#> t2     0     0     1     0
#> t3     0     0     0     1
#> t4    -1    -1    -1    -1
#> 
#> Sum of squares and products for the hypothesis:
#>          time1    time2    time3    time4
#> time1 524.9333 284.3333 507.0667 534.3333
#> time2 284.3333 184.1333 227.4667 396.3333
#> time3 507.0667 227.4667 563.7333 348.6667
#> time4 534.3333 396.3333 348.6667 923.3333
#> 
#> Multivariate Tests: group:time
#>                  Df test stat approx F num Df den Df     Pr(>F)    
#> Pillai            2  1.808841 23.65627      8     20 1.3857e-08 ***
#> Wilks             2  0.008458 22.21451      8     18 7.9942e-08 ***
#> Hotelling-Lawley  2 20.599670 20.59967      8     16 4.8560e-07 ***
#> Roy               2 13.375814 33.43953      4     10 9.2057e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
#> 
#>              Sum Sq num Df Error SS den Df    F value    Pr(>F)    
#> (Intercept) 1155433      1   946.48     12 14649.2234 < 2.2e-16 ***
#> group           912      2   946.48     12     5.7829   0.01743 *  
#> time           2336      4   263.12     48   106.5576 < 2.2e-16 ***
#> group:time      838      8   263.12     48    19.1006 1.621e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> 
#> Mauchly Tests for Sphericity
#> 
#>            Test statistic p-value
#> time              0.29307 0.17766
#> group:time        0.29307 0.17766
#> 
#> 
#> Greenhouse-Geisser and Huynh-Feldt Corrections
#>  for Departure from Sphericity
#> 
#>             GG eps Pr(>F[GG])    
#> time       0.67869  < 2.2e-16 ***
#> group:time 0.67869  4.263e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>              HF eps   Pr(>F[HF])
#> time       0.896371 4.649080e-21
#> group:time 0.896371 2.041727e-11
# 事后检验 不调整与SPSS一样

## 看各时间点处理组成对比较
emmeans::emmeans(SBP_aovez ,  ~ group|time, adjust = "none") %>% pairs()
#> time = t0:
#>  contrast estimate   SE df t.ratio p.value
#>  A - B        -0.2 2.43 12  -0.082  0.9963
#>  A - C        -5.2 2.43 12  -2.137  0.1238
#>  B - C        -5.0 2.43 12  -2.055  0.1415
#> 
#> time = t1:
#>  contrast estimate   SE df t.ratio p.value
#>  A - B        -7.4 3.13 12  -2.364  0.0847
#>  A - C       -10.6 3.13 12  -3.386  0.0139
#>  B - C        -3.2 3.13 12  -1.022  0.5776
#> 
#> time = t2:
#>  contrast estimate   SE df t.ratio p.value
#>  A - B         0.4 2.95 12   0.136  0.9899
#>  A - C        -0.2 2.95 12  -0.068  0.9975
#>  B - C        -0.6 2.95 12  -0.204  0.9774
#> 
#> time = t3:
#>  contrast estimate   SE df t.ratio p.value
#>  A - B        -2.4 3.11 12  -0.771  0.7272
#>  A - C       -16.8 3.11 12  -5.396  0.0004
#>  B - C       -14.4 3.11 12  -4.625  0.0016
#> 
#> time = t4:
#>  contrast estimate   SE df t.ratio p.value
#>  A - B       -14.4 2.50 12  -5.771  0.0002
#>  A - C        -9.8 2.50 12  -3.927  0.0053
#>  B - C         4.6 2.50 12   1.843  0.1975
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates



## 看各处理组时间点成对比较
emmeans::emmeans(SBP_aovez , pairwise ~ time| group, adjust = "none")
#> $emmeans
#> group = A:
#>  time emmean   SE df lower.CL upper.CL
#>  t0      121 1.72 12      117      125
#>  t1      112 2.21 12      108      117
#>  t2      118 2.08 12      114      123
#>  t3      126 2.20 12      121      131
#>  t4      121 1.76 12      117      125
#> 
#> group = B:
#>  time emmean   SE df lower.CL upper.CL
#>  t0      121 1.72 12      117      125
#>  t1      120 2.21 12      115      125
#>  t2      118 2.08 12      113      123
#>  t3      128 2.20 12      123      133
#>  t4      135 1.76 12      131      139
#> 
#> group = C:
#>  time emmean   SE df lower.CL upper.CL
#>  t0      126 1.72 12      122      130
#>  t1      123 2.21 12      118      128
#>  t2      119 2.08 12      114      123
#>  t3      143 2.20 12      138      147
#>  t4      131 1.76 12      127      134
#> 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#> group = A:
#>  contrast estimate    SE df t.ratio p.value
#>  t0 - t1       8.6 1.490 12   5.772  0.0001
#>  t0 - t2       2.6 1.320 12   1.964  0.0732
#>  t0 - t3      -4.8 2.060 12  -2.333  0.0379
#>  t0 - t4       0.2 1.680 12   0.119  0.9074
#>  t1 - t2      -6.0 0.913 12  -6.573  <.0001
#>  t1 - t3     -13.4 1.060 12 -12.624  <.0001
#>  t1 - t4      -8.4 1.530 12  -5.507  0.0001
#>  t2 - t3      -7.4 1.460 12  -5.066  0.0003
#>  t2 - t4      -2.4 1.340 12  -1.789  0.0989
#>  t3 - t4       5.0 1.630 12   3.062  0.0099
#> 
#> group = B:
#>  contrast estimate    SE df t.ratio p.value
#>  t0 - t1       1.4 1.490 12   0.940  0.3659
#>  t0 - t2       3.2 1.320 12   2.417  0.0325
#>  t0 - t3      -7.0 2.060 12  -3.402  0.0052
#>  t0 - t4     -14.0 1.680 12  -8.317  <.0001
#>  t1 - t2       1.8 0.913 12   1.972  0.0721
#>  t1 - t3      -8.4 1.060 12  -7.914  <.0001
#>  t1 - t4     -15.4 1.530 12 -10.096  <.0001
#>  t2 - t3     -10.2 1.460 12  -6.983  <.0001
#>  t2 - t4     -17.2 1.340 12 -12.820  <.0001
#>  t3 - t4      -7.0 1.630 12  -4.287  0.0011
#> 
#> group = C:
#>  contrast estimate    SE df t.ratio p.value
#>  t0 - t1       3.2 1.490 12   2.148  0.0529
#>  t0 - t2       7.6 1.320 12   5.740  0.0001
#>  t0 - t3     -16.4 2.060 12  -7.971  <.0001
#>  t0 - t4      -4.4 1.680 12  -2.614  0.0226
#>  t1 - t2       4.4 0.913 12   4.820  0.0004
#>  t1 - t3     -19.6 1.060 12 -18.465  <.0001
#>  t1 - t4      -7.6 1.530 12  -4.982  0.0003
#>  t2 - t3     -24.0 1.460 12 -16.432  <.0001
#>  t2 - t4     -12.0 1.340 12  -8.944  <.0001
#>  t3 - t4      12.0 1.630 12   7.348  <.0001

8.3.3 分部计算

\[ \begin{aligned} SS_总&= SS_{受试对象间}+SS_{受试对象内} \\ &=（SS_{处理方法}+SS_{个体间差异}）+（SS_{时间}+SS_{处理与时间交互}+SS_{个体内差异}） \end{aligned} \]

\[ \begin{aligned} \nu_总&= \nu_{受试对象间}+\nu_{受试对象内} \\ &=（\nu_{处理方法}+\nu_{个体差异}）+（\nu_{时间}+\nu_{处理与时间交互}+\nu_{个体差异}） \end{aligned} \]

Page 86-87 变异 SS v MS F CP
受试对象间 处理 912.24 2 456.12 5.78 0.0174
个体差异 946.48 12 78.87
1858.72
受试对象内 时间 2336.45 4 584.11 106.56 <0.0001
处理与时间交互 837.63 8 104.70 19.10 <0.0001
个体差异 263.12 48 5.48
3437.2
总 5295.92 3×15-1=74

8.3.3.1 SS总

\[ SS_总=\sum_i^n{（X_i-\bar X）^2} \]

其中，n是观测值的总数，Xi 为每个观测值，\(\bar X\) 是所有观测值的均值。

Show the code

SBP_data <- c(SBP$t0,SBP$t1,SBP$t2,SBP$t3,SBP$t4)

SBP_mean <- mean(SBP_data)

SS_total <- sum((SBP_data-SBP_mean)^2)

SS_total    
#> [1] 5295.92

8.3.3.2 SS受试对象间

\[ SS_{受试对象间}=\sum_{j=1}^m n_{.j}(\bar X_{.j}-\bar X)^2 \]

其中，m是受试者数量（15），n.j 是第j 个受试对象的观测值数量，\(\bar X.j\)是第j 个受试对象的观测值的均值，\(\bar X\) 是所有观测值的均值。

Show the code
id_mean <- SBP |> dplyr::select(c(-1,-2)) |> rowMeans()

SS_between <- 0
for (i in 1:nrow(SBP)) {
    SS_between <- SS_between + 5*(id_mean[i]-SBP_mean)^2
}
SS_between
#> [1] 1858.72
Show the code
group__summary <- SBP_long |> group_by(group) |> 
    summarise(n=n(),mean=mean(SBP),sum=sum(SBP))
group__summary 
#> # A tibble: 3 × 4
#>   group     n  mean   sum
#>   <fct> <int> <dbl> <dbl>
#> 1 A        25  120.  2992
#> 2 B        25  124.  3112
#> 3 C        25  128.  3205
SS处理

\[ SS_{处理}= \sum_i^k n_k (\bar X_k - \bar X)^2 \]

其中，k为不同处理组的组数，nk 为第 k 处理组观测值的总数，\(\bar X_k\) 为第k 处理组观测值的均数，\(\bar X\) 是所有观测值的均值。

Show the code
SS_处理 <- 25*( (group__summary$mean[1] - SBP_mean )^2 + 
                       (group__summary$mean[2] - SBP_mean )^2 +
                       (group__summary$mean[3] - SBP_mean )^2) 
SS_处理
#> [1] 912.24
Show the code
SS_between_error <- SS_between-SS_处理

SS_between_error
#> [1] 946.48

8.3.3.3 SS受试对象内

\[ SS_{受试对象内}=SS_总-SS_{受试对象间} \]

Show the code
SS_within <- SS_total-SS_between
SS_within
#> [1] 3437.2
SS时间

\[ SS_{时间}=\sum _{t=1}^T n_t （\bar X_{.t}-\bar X）^2 \]

其中，T是时间点的数量，nt 是在时间点t的观测值数量，\bar X.t 是在时间点 t 的均值。

Show the code
t_mean <- SBP |> dplyr::select(c(-1,-2)) |> colMeans()

SS_time <- 0
for (i in seq_along(t_mean)) {
    
    SS_time <- SS_time + 15*(t_mean[i]-SBP_mean)^2
    
}
names(SS_time) <- "SS_time"
SS_time
#>  SS_time 
#> 2336.453
SS处理与时间交互

\[ SS_{处理与时间交互}=\sum_{i=1}^{k}\sum_{j=1}^{T}n_{ij}\left (\bar X_{ij.} -\bar X_{i..}-\bar X_{.j.} + \bar X_{...} \right )^2 \]

其中，

  • k 是处理方法的数量，3

  • T 是时间点的数量，5

  • nij 是第 i 个处理方法和第 j 个时间点的观测次数，25/5=5

  • Xij . 是第 i 个处理方法、第 j 个时间点的观测值的平均值

  • Xi . . 是第 i 个处理方法观测值的平均值，不考虑时间点

  • X. j . 是第 j 个时间点观测值的平均值，不考虑处理方法

  • X. . . 是所有观测值的平均值，不考虑处理方法和时间点

Show the code
interaction_effect_summary <- SBP_long |>
    summarise(
        #n = n(),
        mean = mean(SBP),
       # sum = sum(SBP),
        .by = c(group , time)
    )
interaction_effect_mean <- interaction_effect_summary |> pivot_wider(names_from = c("time"),values_from = "mean")



SS_interaction_effect <- 0
for (i in 1:3) {
    for (j in 1:5) {
        SS_interaction_effect <- SS_interaction_effect + 5 * (interaction_effect_mean[i, j + 1] -group__summary$mean[i] - t_mean[j] + SBP_mean) ^ 2
    }
    
}
names(SS_interaction_effect) <- "SS_interaction_effect"
SS_interaction_effect
#>   SS_interaction_effect
#> 1              837.6267
Show the code
SS_within_error <- SS_within-SS_time-SS_interaction_effect

names(SS_within_error) <- "SS_within_error"
SS_within_error 
#>   SS_within_error
#> 1          263.12

8.4 球形度检验 sphericity

重复测量方差分析假设相关条件（或组水平）的所有组合之间的差异方差相等。这被称为球形度假设。

R中球形度检验

球度仅针对具有两个以上水平的变量进行评估，因为球形度必然适用于只有两个水平的条件。

违反球形度假设可能会扭曲方差计算，从而导致更宽松的重复测量方差分析检验（即 I 类错误率增加）。在这种情况下，必须根据违反球形度的程度适当校正重复测量方差分析。文献中使用了两种常见的校正：Greenhouse-Geisser epsilon （GGe） 和 Huynh-Feldt epsilon （HFe）。

球形度的 Mauchly 检验用于评估是否满足球形度的假设。使用rstatix::anova_test() 时，会自动报告此问题。尽管该方法受到严厉批评，通常无法检测到小样本中的球形偏离，而在大样本中则过度检测到它们，但它仍然是一种常用的方法。

8.4.1 计算步骤

具体操作步骤如下：

  1. 计算相关组的每个组合之间的差异

  2. 计算每个组差的方差

8.4.2 案例

8.4.2.1 rstatix::anova_test()

Show the code

aov_SBP <- rstatix::anova_test(data = SBP_long,
           dv = SBP,
           wid = id,
           within = time,
           between = group
           )

aov_SBP
#> ANOVA Table (type II tests)
#> 
#> $ANOVA
#>       Effect DFn DFd       F        p p<.05   ges
#> 1      group   2  12   5.783 1.70e-02     * 0.430
#> 2       time   4  48 106.558 3.02e-23     * 0.659
#> 3 group:time   8  48  19.101 1.62e-12     * 0.409
#> 
#> $`Mauchly's Test for Sphericity`
#>       Effect     W     p p<.05
#> 1       time 0.293 0.178      
#> 2 group:time 0.293 0.178      
#> 
#> $`Sphericity Corrections`
#>       Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
#> 1       time 0.679 2.71, 32.58 1.87e-16         * 0.896 3.59, 43.03 4.65e-21
#> 2 group:time 0.679 5.43, 32.58 4.26e-09         * 0.896 7.17, 43.03 2.04e-11
#>   p[HF]<.05
#> 1         *
#> 2         *

输出是一个包含三个表的列表：

  1. 方差分析结果显示标有（generalized eta squared，ges）的列上的 p 值和效应大小;效应大小本质上是由于受试者内因素而忽略受试者效应而导致的变异量。

  2. Mauchly Sphericity检验。仅报告具有 >2 水平的变量或效应，因为球形度必然适用于只有 2 个水平的效应。原假设是组差的方差相等。因此，显著的 p 值 （p < 0.05） 表示组差的方差不相等。

  3. 球度校正结果，以防我们无法维持球形度假设。提供了文献中使用的两种常见校正：Greenhouse-Geisser epsilon （GGe） 和 Huynh-Feldt epsilon （HFe） 及其相应的 p 值

8.4.2.2 ez::ezANOVA()

Show the code
library(ez)

# 进行重复测量方差分析
aov_ez <- ezANOVA(data = SBP_long,
                   dv = SBP,
                   wid = id,
                   within = time,
                   between = group,
                   detailed = TRUE)

# 查看球形检验结果
print(aov_ez)
#> $ANOVA
#>        Effect DFn DFd          SSn    SSd            F            p p<.05
#> 1 (Intercept)   1  12 1155433.0800 946.48 14649.223396 6.784700e-20     *
#> 2       group   2  12     912.2400 946.48     5.782943 1.743351e-02     *
#> 3        time   4  48    2336.4533 263.12   106.557616 3.017101e-23     *
#> 4  group:time   8  48     837.6267 263.12    19.100638 1.621310e-12     *
#>         ges
#> 1 0.9989542
#> 2 0.4299287
#> 3 0.6588884
#> 4 0.4091519
#> 
#> $`Mauchly's Test for Sphericity`
#>       Effect         W         p p<.05
#> 3       time 0.2930746 0.1776608      
#> 4 group:time 0.2930746 0.1776608      
#> 
#> $`Sphericity Corrections`
#>       Effect      GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
#> 3       time 0.678693 1.867336e-16         * 0.896371 4.649080e-21         *
#> 4 group:time 0.678693 4.262529e-09         * 0.896371 2.041727e-11         *

8.4.3 当满足球形度假设时

球形度的 Mauchly 检验不显著 （p > 0.05）;这表明，受试者内因素水平之间的差异方差是相等的。因此，我们可以假设协方差矩阵的球形度，并解释方差分析表中可用的标准输出。

Show the code
# Display ANOVA table
aov_SBP$ANOVA
#>       Effect DFn DFd       F        p p<.05   ges
#> 1      group   2  12   5.783 1.70e-02     * 0.430
#> 2       time   4  48 106.558 3.02e-23     * 0.659
#> 3 group:time   8  48  19.101 1.62e-12     * 0.409

  • F表示我们正在与 F 分布（F 检验）进行比较; 分别表示 time 和 Error（time） 的自由度; 表示得到的 F 统计量值(2, 18)81.8

  • p指定 p 值

  • ges（广义 eta 平方，etaov_SBP[g]）是效应大小（由于受试者内因素引起的变异量）

8.4.4 当违反球形度假设时

如果数据违反了球形度假设（即 Mauchly 检验，p <= 0.05），则应解释表sphericity corrections中的结果，其中对自由度进行了调整，这会影响检验的统计显著性（即 p 值）。校正通过乘法和校正估计值（Greenhouse-Geisser （GG） 和 Huynh-Feldt （HF） ε 值）来应用。

Note

epsilon 提供了球形度被侵犯的程度的度量。值为 1 表示不偏离球形度（组差的所有方差均相等）。违反球形度会导致 epsilon 值低于 1。epsilon 离 1 越远，违规越严重。

可以看出，即使在球形度校正（p[GG] < 0.001，p[HF] < 0.001）之后，平均自尊得分在不同时间点仍存在统计学差异。

在两种球形度校正方法中，Huynh-Feldt 校正被认为是最不保守的（高估了 epsilon），而 Greenhouse-Geisser 被认为是更保守的（当 epsilon 接近 1 时低估了 epsilon）。

一般建议使用 Greenhouse-Geisser 校正，特别是当 epsilon < 0.75 时。在 epsilon 大于 0.75 的情况下，一些统计学家建议使用 Huynh-Feldt 校正 （Girden 1992）。