# 重复测量方差分析 - **总体均值**： 代表所有观测值的平均值，是固定效应的一部分。- **组间效应**： 反映不同组（如治疗组和对照组）之间的差异，是固定效应的一部分。- **时间效应(**： 反映不同时间点或条件（如测量前、测量中、测量后）之间的差异，是固定效应的一部分。- **组间和时间效应的交互作用**： 反映不同组在不同时间点或条件下的差异，是固定效应的一部分。- **受试者效应**： 反映个体差异，考虑每个受试者的独特性，通常作为随机效应处理。- **时间效应与受试者效应的交互作用**： 反映个体在不同时间点或条件下的反应差异，通常作为随机效应处理。- **误差项**： 表示观测值中的随机误差或未解释的变异。## 假设条件 - 独立性假设，即不同受试者之间的观测值是独立的。- 观测值在每个时间点或条件下的正态分布。- 各组内的方差齐性（sphericity），即组内相关性结构是已知的，通常假设为复合对称（compound symmetry）。**协方差矩阵的球形检验（W=1）**- **复合对称假设**：RM-ANOVA假设误差项具有复合对称性，即组内相关性和方差齐性。如果复合对称假设不成立，结果可能不准确，可以使用更为灵活的线性混合模型（LMM）或广义估计方程（GEE）进行分析。- **检验方法**：RM-ANOVA的主要检验方法包括F检验，用于检验组效应、时间效应和交互作用是否显著。## 单因素重复测量方差分析 ### 原理公式 ```{r} <- 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 ,<- |> pivot_longer (cols = - 1 ,names_to = "week" ,values_to = "ALT" ) |> mutate (week = factor (week))# 假设检验 shapiro.test (df$ zero)shapiro.test (df$ ten)shapiro.test (df$ twenty)shapiro.test (df$ thirty)bartlett.test (df[- 1 ])# 总变异 <- sum (df_long$ ALT)<- sum (df_long$ ALT^ 2 )<- xij_sum^ 2 / 40 <- xij_square_sum- C# 受试者间 <- rowSums (df[,- 1 ])<- sum (xi._sum^ 2 )/ 4 - C #行和平方和 <- SS_B/ 9 # 不同时间点 <- colSums (df[,- 1 ])<- sum (x.j_sum^ 2 )/ 10 - C #列和平方和 <- SS_W/ 3 <- SS_T- SS_B- SS_W<- SS_E/ (9 * 3 )<- MS_W/ MS_E``` ### 线性混合模型`nlme::lme()` ```{r} library (nlme)<- lme (fixed= ALT ~ week, random = ~ 1 | id/ week, data = df_long)summary (model)<- anova (model)``` ### 事后检验 ```{r} library (multcomp)<- glht (model, linfct = mcp (week = c ("ten - zero = 0 " ,"twenty - zero == 0" ,"thirty - zero = 0 " )))summary (glht_result)``` ### 第二个示例 ```{r} <- tribble (~ Participant,~ A1,~ A2,~ A3,"P1" ,5 ,6 ,7 ,"P2" ,7 ,13 ,13 ,"P3" ,2 ,4 ,6 ,"P4" ,6 ,9 ,12 ``` ```{r} <- 3-1 <- 4-1 <- df_A * df_P<- df_A + df_P + df_error# 或者 <- 3 * 4-1 ``` ## 含组间因子和组内因子的重复测量 \\ <https://personality-project.org/r/r.guide/r.anova.html#oneway> ### 第一个案例 ```{r} <- read_rds ("data/SLAMF6+-PBS.rds" )%>% :: datatable ()``` ```{r} #| code-fold: show # 2x2 mixed: 独立变量（被试间） : age # 独立变量（被试内） : time 依赖变量: value # aov_age_time <- aov(value ~ age * time + Error(subject/time), # data = data_long) # summary(aov_age_time) # 两个被试内变量 #aov.ww <- aov(y ~ w1*w2 +Error(subject/(w1*w2)), data=data_long) # 1个被试间变量，两个被试内变量 #aov.bww <- aov(y ~b1*w1*w2 + Error(subject/(w1*w2)) + b1,data=data_long) # 两个被试间变量，一个被试内变量 # aov.bww <- aov(y ~ b1*b2*w1 + Error(subject/(w1)) + b1*b2, data=data_long) ``` #### `stats::aov()` ```{r} # 第一种 有混合效应无法事后检验 <- aov (volume ~ method* Days+ Error (id/ Days),data = df)summary (aov_1)``` #### `nlme::lme()` ```{r} # 第二种 最好 library (nlme)<- lme (volume ~ method* Days, random = ~ 1 | id/ Days, data = df)summary (model)<- anova (model)#成对比较 library (emmeans)<- emmeans (model, ~ method) print (method_means) <- pairs (method_means) print (method_comparisons)``` #### `afex::aov_*()` ```{r} library ("afex" ) # needed for ANOVA functions. library ("emmeans" ) # emmeans must now be loaded explicitly for follow-up tests. library ("multcomp" ) # for advanced control for multiple testing/Type 1 errors. library ("ggplot2" ) # for customizing plots. <- aov_ez (id = "id" ,dv = "volume" , df, between = "method" , within = c ("Days" ))aov_car (volume ~ method + Error (id/ Days), data = df)aov_4 (volume ~ method + (Days| id), data = df)# 事后检验 <- emmeans (a1, ~ method)pairs (m1)``` ### 第二个案例 ```{r} # 医学统计学 人卫版 第7版 <- read_delim ("data/麻醉诱导时相.txt" )<- df |> pivot_longer (cols = starts_with ("t" ),names_to = "time" ,values_to = "BP" ) |> mutate (id= factor (id),group= factor (group),time= factor (time)%>% :: datatable ()``` #### `nlme::lme()` ```{r} # 第一种方法 library (nlme)<- lme (BP ~ group* time, random = ~ 1 | id/ time, data = df_long)summary (model)<- anova (model)#成对比较 library (emmeans)<- emmeans (model, ~ group) print (method_means) <- pairs (method_means) print (method_comparisons)``` #### `rstatix::anova_test()` ```{r} library (rstatix)<- anova_test (data = df_long,dv = BP,wid = id,within = time,between = group``` #### `afex::aov_*()` <https://cran.r-project.org/web/packages/afex/vignettes/afex_anova_example.html> ```{r} library ("afex" ) # needed for ANOVA functions. library ("emmeans" ) # emmeans must now be loaded explicitly for follow-up tests. library ("multcomp" ) # for advanced control for multiple testing/Type 1 errors. library ("ggplot2" ) # for customizing plots. <- aov_ez ("id" , "BP" , df_long, between = "group" , within = c ("time" ))aov_car (BP ~ group + Error (id/ time), data = df_long)aov_4 (BP ~ group + (time| id), data = df_long)# 事后检验 <- emmeans (a3, ~ group)pairs (m3)summary (as.glht (pairs (m1)), test= adjusted ("fdr" ))<- afex_plot (a3, x = "time" ,# trace = "BP", panel = "group" , error = "none" ,mapping = c ("color" , "fill" ),data_geom = geom_boxplot, data_arg = list (width = 0.4 ),point_arg = list (size = 1.5 ), line_arg = list (size = 1 ))``` \\ \\ \< 0.0001 |\< 0.0001 |#### SS~总~ ```{r} options (digits = 4 )<- c (df$ t0,df$ t1,df$ t2,df$ t3,df$ t4)<- mean (BP) <- sum ((BP- BP_mean)^ 2 )``` #### SS~受试对象间~ \\ bar X~.j~ 是第j 个受试对象的观测值的均值，\\ bar X 是所有观测值的均值。```{r} <- df |> dplyr:: select (c (- 1 ,- 2 )) |> rowMeans ()<- 0 for (i in 1 : nrow (df)) {<- SS_between + 5 * (id_mean[i]- BP_mean)^ 2 ``` ```{r} <- df_long |> group_by (group) |> summarise (n= n (),mean= mean (BP),sum= sum (BP))``` ##### SS~处理~ ```{r} <- 25 * ( (group__summary$ mean[1 ] - BP_mean )^ 2 + $ mean[2 ] - BP_mean )^ 2 + $ mean[3 ] - BP_mean )^ 2 ) ``` ```{r} <- SS_between- SS_处理``` #### SS~受试对象内~ ```{r} <- SS_total- SS_between``` ##### SS~时间~ \\ bar X~.t~ 是在时间点 t 的均值。```{r} <- df |> dplyr:: select (c (- 1 ,- 2 )) |> colMeans ()<- 0 for (i in seq_along (t_mean)) {<- SS_time + 15 * (t_mean[i]- BP_mean)^ 2 names (SS_time) <- "SS_time" ``` ##### SS~处理与时间交互~ - k 是处理方法的数量，3- T 是时间点的数量，5- n~ij~ 是第 i 个处理方法和第 j 个时间点的观测次数，25/5=5- \bar Xij ~.~ 是第 i 个处理方法、第 j 个时间点的观测值的平均值- \bar X~i . .~ 是第 i 个处理方法观测值的平均值，不考虑时间点- \bar X~. j .~ 是第 j 个时间点观测值的平均值，不考虑处理方法- \bar X~. . .~ 是所有观测值的平均值，不考虑处理方法和时间点```{r} <- df_long |> summarise (#n = n(), mean = mean (BP),# sum = sum(BP), .by = c (group , time)<- interaction_effect_summary |> pivot_wider (names_from = c ("time" ),values_from = "mean" )<- 0 for (i in 1 : 3 ) {for (j in 1 : 5 ) {<- SS_interaction_effect + 5 * (interaction_effect_mean[i, j + 1 ] - group__summary$ mean[i] - t_mean[j] + BP_mean) ^ 2 names (SS_interaction_effect) <- "SS_interaction_effect" ``` ```{r} <- SS_within- SS_time- SS_interaction_effectnames (SS_within_error) <- "SS_within_error" ``` ## 球形度检验 sphericity [ R中球形度检验 ](https://www.datanovia.com/en/lessons/mauchlys-test-of-sphericity-in-r/) `rstatix::anova_test()` 时，会自动报告此问题。尽管该方法受到严厉批评，通常无法检测到小样本中的球形偏离，而在大样本中则过度检测到它们，但它仍然是一种常用的方法。```{r} library (tidyverse)library (ggpubr)library (rstatix)``` ```{} ``` ### 计算 sphericity 和 Mauchly 检验 1. 计算相关组的每个组合之间的差异2. 计算每个组差的方差 ```{r} <- df |> mutate (` t0-t1 ` = t0- t1,` t0-t2 ` = t0- t2,` t0-t3 ` = t0- t3,` t0-t4 ` = t0- t4,` t1-t2 ` = t1- t2,` t1-t3 ` = t1- t3,` t1-t4 ` = t1- t4,` t2-t3 ` = t2- t3,` t2-t4 ` = t2- t4,` t3-t4 ` = t3- t4,|> select (8 : 17 )|> map (var)### 单因素 <- df |> select (3 : 7 ) |> as.matrix ()mauchly.test (lm (t~ 1 ),X = ~ 1 )### 双因素 <- df |> group_split (group)1 : 5 ,3 : 7 ]<- as.matrix (cbind (df_split[[1 ]][3 : 7 ],df_split[[2 ]][3 : 7 ],df_split[[3 ]][3 : 7 ]))<- ordered (rep (1 : 5 ,3 ))<- factor (rep (c ("A" ,"B" ,"C" ),each= 5 ))mauchly.test (lm (df2~ 1 ),M= ~ group+ times,X= ~ times)``` ### 方差分析 #### `rstatix::anova_test()` ```{r} <- rstatix:: anova_test (data = df_long,dv = BP,wid = id,within = time,between = group``` 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 值#### `ez::ezANOVA()` ```{r} library (ez)# 进行重复测量方差分析 <- ezANOVA (data = df_long,dv = BP,wid = id,within = time,between = group,detailed = TRUE )# 查看球形检验结果 print (aov_ez)``` ### **当满足球形度假设时** \> 0.05）;这表明，受试者内因素水平之间的差异方差是相等的。因此，我们可以假设协方差矩阵的球形度，并解释方差分析表中可用的标准输出。```{r} # Display ANOVA table $ ANOVA``` - `F` 表示我们正在与 F 分布（F 检验）进行比较; 分别表示 time 和 Error（time） 的自由度; 表示得到的 F 统计量值`(2, 18)81.8` - `p` 指定 p 值- `ges` （广义 eta 平方，eta2\[ g\] ）是效应大小（由于受试者内因素引起的变异量）### **当违反球形度假设时** \< = 0.05），则应解释表`sphericity corrections` 中的结果，其中对自由度进行了调整，这会影响检验的统计显著性（即 p 值）。校正通过乘法和校正估计值（Greenhouse-Geisser （GG） 和 Huynh-Feldt （HF） ε 值）来应用。```{r} $ ` Sphericity Corrections ` ``` \[ GG\] \< 0.001，p\[ HF\] \< 0.001）之后，平均自尊得分在不同时间点仍存在统计学差异。\< 0.75 时。在 epsilon 大于 0.75 的情况下，一些统计学家建议使用 Huynh-Feldt 校正 （Girden 1992）。```{r} # correction = "auto" get_anova_table (a2)``` ```{r} # correction = "GG" get_anova_table (a2, correction = "GG" )``` ### post-hoc test ```{r} |> pairwise_t_test (BP~ time,paired = T,p.adjust.method = "bonferroni" )``` 
