如何在R中的多列中执行三组、三个时间点和一个因变量的双向混合方差分析
How to perform a two-way mixed ANOVA with three groups, three time points and one dependent variable in multiple columns in R
我从大脑中的几个感兴趣区域 (fMRI) 中提取了时间序列，并且在对应于大脑中两个节点之间相关性的列下为每个受试者添加了成对相关 (Fisher-Z) 值（对于示例：stim_lvis3，stim = 刺激部位和 lvis3= 左视觉网络 3)。现在，我想对该数据集执行方差分析以查看效果和 between/within 组差异（3 组 x 3 个时间点）。我的数据已经是长格式了。
*groups= ctbs [10 个科目 x 3 个时间点]，itbs = [10 个科目 x 3 个时间点]，以及假 [10 个科目 x 3 个时间点]
鉴于我有 12 列具有连接值 (stim_lvis3...stim_rpcc1)，关于如何完成此操作的任何建议。例如，我无法按时间和组对数据进行箱线图？
如何在这种情况下对每个组在特定时间点的所有 12 列执行双向混合方差分析，然后在每个时间点比较组？
我将主题、时间和组别转换为因素
tbs %>%
group_by(time, group) %>%
get_summary_stats(stim_lVis3, type = "mean_sd")
Error in tbs(.) : could not find function "tbs"
bxp <- ggboxplot(
tbs, x = "time", y = "stim_lvis3",
color = "group", palette = "jco"
)
bxp
Error in FUN(X[[i]], ...) : object 'stim_lvis3' not found
欢迎来到 SO！看起来你是新手，但将来要快速获得很好的答案，请确保你以可用的格式包含样本数据（即 dput(head(myData)) 的输出）。检查一下：making R reproducible questions.
我知道有两种方法可以在方差分析内和方差分析之间完成。更容易实现的版本来自包ez。包 jmv 提供了一个更复杂的 write-up 但你也有很大的控制权。我确定 ez 的版本还有更多内容，但我并没有太多地使用那个包。
我创建了一些数据来模拟您正在使用的内容。
library(tidyverse)
library(jmv)
library(ez)
set.seed(35)
df1 <- data.frame(subject = factor(rep(1:15, each = 3)),
time = factor(rep(c(1:3), 15)),
group = factor(rep(1:3, each = 15)),
stim_IVis3 = rnorm(45, .5, .15),
stim_IVis4 = rnorm(45, .51, .125))
head(df1)
summary(df1)
要使用 ez，这是一个非常简单的实现。虽然，我找不到允许多个因变量的选项。本质上，您要么需要 pivot_long，要么可以使用 jmv。
对于这种方法，您不会得到 post 临时比较；效应大小是广义的 η2.
ezANOVA(df1, dv = stim_IVis3, wid = subject, within = time,
between = group, detailed = T)
# $ANOVA
# Effect DFn DFd SSn SSd F p p<.05 ges
# 1 (Intercept) 1 12 12.58700023 0.3872673 390.0252306 1.616255e-10 * 0.92579702
# 2 group 2 12 0.05853169 0.3872673 0.9068417 4.297644e-01 0.05483656
# 3 time 2 24 0.05417372 0.6215855 1.0458491 3.668654e-01 0.05096178
# 4 group:time 4 24 0.06774352 0.6215855 0.6539102 6.298074e-01 0.06292379
#
# $`Mauchly's Test for Sphericity`
# Effect W p p<.05
# 3 time 0.9914977 0.9541232
# 4 group:time 0.9914977 0.9541232
#
# $`Sphericity Corrections`
# Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
# 3 time 0.9915694 0.3664558 1.187039 0.3668654
# 4 group:time 0.9915694 0.6286444 1.187039 0.6298074
现在要使用 jmv，您需要将数据旋转得更宽。您的 within-subject 数据需要在单独的列中。由于 time 是用 stim... 列中的 within-subject 值表示的因子，因此这就是您需要调整的内容。
df2 <- df1 %>%
pivot_wider(id_cols = c(subject, group),
names_from = time,
values_from = starts_with("stim"),
names_vary = "fastest")
head(df2)
# # A tibble: 6 × 8
# subject group stim_IVis3_1 stim_IVis3_2 stim_IVis3_3 stim_IVis4_1 stim_IVis4_2 stim_IVis4_3
# <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 1 1 0.497 0.505 0.601 0.417 0.776 0.660
# 2 2 1 0.584 0.408 0.737 0.454 0.670 0.650
# 3 3 1 0.741 0.399 0.450 0.306 0.665 0.577
# 4 4 1 0.223 0.450 0.350 0.522 0.342 0.417
# 5 5 1 0.717 0.661 0.581 0.601 0.700 0.686
# 6 6 2 0.534 0.780 0.443 0.420 0.565 0.739
现在您已准备好 jmv。 bs 等同于 between-subjects； rm 等于重复测量。
fit = anovaRM(data = df2,
ss = "3", # type of SS (1, 2, or 3)
bs = list("group"), # between subjects (links with other parameters this way)
bsTerms = list("group"),
rm = list(list(label = "stim_IVs", # within subjects
levels = c(names(df2)[3:8]))), # could also write c("stim_IVs3_1", "stim_IVs3_2", "stim_IVs3_3")
rmCells = list(list(measure = names(df2)[3], # could also write "stim_IVs3_1"
cell = names(df2)[3]), # the group associated in 'rm' level (I used the column name)
list(measure = names(df2)[4],
cell = names(df2)[4]),
list(measure = names(df2)[5],
cell = names(df2)[5]),
list(measure = names(df2)[6],
cell = names(df2)[6]),
list(measure = names(df2)[7],
cell = names(df2)[7]),
list(measure = names(df2)[8],
cell = names(df2)[8])
),
rmTerms = list("stim_IVs"), # groups variable for repeated/within measures
emMeans = list(list("stim_IVs", "group")), # all grouping vars in the ANOVA (for the em tables)
emmPlots = T, # show em plots
emmTables = T, # show em tables
effectSize = c("omega", "partEta"), # multiple options here, see help
spherTests = T, # use correction
spherCorr = c("none", "GG", "HF"), # multiple options here, see help
leveneTest = T, # check homogeneity (p GREATER than .05 is good)
qq = T, # plot normality validation qq plot
postHoc = list("group",c("group","stim_IVs"),"stim_IVs"),
postHocCorr = "tukey") # use TukeyHSD
如果您打印 fit，您将获得 吨 的信息。除了 ezANOVA 提供的内容外，还提供了 post 每个组、时间、因变量以及每个组的混合的临时比较；统计假设的结果：Levene 的同质性和 QQ 正态残差图；以及估计的边际均值 table 和这些均值的图。
我发现您的数据中的字段比我这里的多。我相信我已经为您提供了足够的基础。
我建议您确定要从数据中学习什么，然后从中选择算法。
如果您有任何问题，请告诉我。
我从大脑中的几个感兴趣区域 (fMRI) 中提取了时间序列，并且在对应于大脑中两个节点之间相关性的列下为每个受试者添加了成对相关 (Fisher-Z) 值（对于示例：stim_lvis3，stim = 刺激部位和 lvis3= 左视觉网络 3)。现在，我想对该数据集执行方差分析以查看效果和 between/within 组差异（3 组 x 3 个时间点）。我的数据已经是长格式了。
*groups= ctbs [10 个科目 x 3 个时间点]，itbs = [10 个科目 x 3 个时间点]，以及假 [10 个科目 x 3 个时间点]
鉴于我有 12 列具有连接值 (stim_lvis3...stim_rpcc1)，关于如何完成此操作的任何建议。例如，我无法按时间和组对数据进行箱线图？
如何在这种情况下对每个组在特定时间点的所有 12 列执行双向混合方差分析，然后在每个时间点比较组？
我将主题、时间和组别转换为因素
tbs %>%
group_by(time, group) %>%
get_summary_stats(stim_lVis3, type = "mean_sd")
Error in tbs(.) : could not find function "tbs"
bxp <- ggboxplot(
tbs, x = "time", y = "stim_lvis3",
color = "group", palette = "jco"
)
bxp
Error in FUN(X[[i]], ...) : object 'stim_lvis3' not found
欢迎来到 SO！看起来你是新手，但将来要快速获得很好的答案，请确保你以可用的格式包含样本数据（即 dput(head(myData)) 的输出）。检查一下：making R reproducible questions.
我知道有两种方法可以在方差分析内和方差分析之间完成。更容易实现的版本来自包ez。包 jmv 提供了一个更复杂的 write-up 但你也有很大的控制权。我确定 ez 的版本还有更多内容，但我并没有太多地使用那个包。
我创建了一些数据来模拟您正在使用的内容。
library(tidyverse)
library(jmv)
library(ez)
set.seed(35)
df1 <- data.frame(subject = factor(rep(1:15, each = 3)),
time = factor(rep(c(1:3), 15)),
group = factor(rep(1:3, each = 15)),
stim_IVis3 = rnorm(45, .5, .15),
stim_IVis4 = rnorm(45, .51, .125))
head(df1)
summary(df1)
要使用 ez，这是一个非常简单的实现。虽然，我找不到允许多个因变量的选项。本质上，您要么需要 pivot_long，要么可以使用 jmv。
对于这种方法，您不会得到 post 临时比较；效应大小是广义的 η2.
ezANOVA(df1, dv = stim_IVis3, wid = subject, within = time,
between = group, detailed = T)
# $ANOVA
# Effect DFn DFd SSn SSd F p p<.05 ges
# 1 (Intercept) 1 12 12.58700023 0.3872673 390.0252306 1.616255e-10 * 0.92579702
# 2 group 2 12 0.05853169 0.3872673 0.9068417 4.297644e-01 0.05483656
# 3 time 2 24 0.05417372 0.6215855 1.0458491 3.668654e-01 0.05096178
# 4 group:time 4 24 0.06774352 0.6215855 0.6539102 6.298074e-01 0.06292379
#
# $`Mauchly's Test for Sphericity`
# Effect W p p<.05
# 3 time 0.9914977 0.9541232
# 4 group:time 0.9914977 0.9541232
#
# $`Sphericity Corrections`
# Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
# 3 time 0.9915694 0.3664558 1.187039 0.3668654
# 4 group:time 0.9915694 0.6286444 1.187039 0.6298074
现在要使用 jmv，您需要将数据旋转得更宽。您的 within-subject 数据需要在单独的列中。由于 time 是用 stim... 列中的 within-subject 值表示的因子，因此这就是您需要调整的内容。
df2 <- df1 %>%
pivot_wider(id_cols = c(subject, group),
names_from = time,
values_from = starts_with("stim"),
names_vary = "fastest")
head(df2)
# # A tibble: 6 × 8
# subject group stim_IVis3_1 stim_IVis3_2 stim_IVis3_3 stim_IVis4_1 stim_IVis4_2 stim_IVis4_3
# <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 1 1 0.497 0.505 0.601 0.417 0.776 0.660
# 2 2 1 0.584 0.408 0.737 0.454 0.670 0.650
# 3 3 1 0.741 0.399 0.450 0.306 0.665 0.577
# 4 4 1 0.223 0.450 0.350 0.522 0.342 0.417
# 5 5 1 0.717 0.661 0.581 0.601 0.700 0.686
# 6 6 2 0.534 0.780 0.443 0.420 0.565 0.739
现在您已准备好 jmv。 bs 等同于 between-subjects； rm 等于重复测量。
fit = anovaRM(data = df2,
ss = "3", # type of SS (1, 2, or 3)
bs = list("group"), # between subjects (links with other parameters this way)
bsTerms = list("group"),
rm = list(list(label = "stim_IVs", # within subjects
levels = c(names(df2)[3:8]))), # could also write c("stim_IVs3_1", "stim_IVs3_2", "stim_IVs3_3")
rmCells = list(list(measure = names(df2)[3], # could also write "stim_IVs3_1"
cell = names(df2)[3]), # the group associated in 'rm' level (I used the column name)
list(measure = names(df2)[4],
cell = names(df2)[4]),
list(measure = names(df2)[5],
cell = names(df2)[5]),
list(measure = names(df2)[6],
cell = names(df2)[6]),
list(measure = names(df2)[7],
cell = names(df2)[7]),
list(measure = names(df2)[8],
cell = names(df2)[8])
),
rmTerms = list("stim_IVs"), # groups variable for repeated/within measures
emMeans = list(list("stim_IVs", "group")), # all grouping vars in the ANOVA (for the em tables)
emmPlots = T, # show em plots
emmTables = T, # show em tables
effectSize = c("omega", "partEta"), # multiple options here, see help
spherTests = T, # use correction
spherCorr = c("none", "GG", "HF"), # multiple options here, see help
leveneTest = T, # check homogeneity (p GREATER than .05 is good)
qq = T, # plot normality validation qq plot
postHoc = list("group",c("group","stim_IVs"),"stim_IVs"),
postHocCorr = "tukey") # use TukeyHSD
如果您打印 fit，您将获得 吨 的信息。除了 ezANOVA 提供的内容外，还提供了 post 每个组、时间、因变量以及每个组的混合的临时比较；统计假设的结果：Levene 的同质性和 QQ 正态残差图；以及估计的边际均值 table 和这些均值的图。
我发现您的数据中的字段比我这里的多。我相信我已经为您提供了足够的基础。
我建议您确定要从数据中学习什么，然后从中选择算法。
如果您有任何问题，请告诉我。