如何在 R 中的多列中执行具有三组、三个时间点和一个因变量的双向混合方差分析

发布于 2025-01-16 23:05:56 字数 925 浏览 1 评论 0原文

我从大脑中的几个感兴趣区域（fMRI）中提取了时间序列，并且在与大脑中两个节点之间的相关性相对应的列下为每个受试者添加了成对相关性（Fisher-Z）值（例如：stim_lvis3、stim = 刺激部位，lvis3= 左视觉网络 3)。现在，我想对此数据集执行方差分析，以查看效果以及组间/组内差异（3 组 x 3 个时间点）。我的数据已经是长格式了。

*groups= ctbs [10 个受试者 x 3 个时间点]、itbs = [10 个受试者 x 3 个时间点] 和 sham [10 个受试者 x 3 个时间点]

考虑到我有 12 列，关于如何做到这一点的任何建议连接值（stim_lvis3...stim_rpcc1）。例如，我无法绘制按时间和组划分的数据箱线图？
在这种情况下，如何在特定时间点对每组的所有 12 列执行双向混合方差分析，然后在每个时间点比较各组？

data in long format

我将主题、时间和组转换为因素

tbs %>%
  group_by(time, group) %>%
  get_summary_stats(stim_lVis3, type = "mean_sd")

tbs(.) 中的错误：找不到函数“tbs”

bxp <- ggboxplot(
  tbs, x = "time", y = "stim_lvis3",
  color = "group", palette = "jco"
)

bxp

FUN(X[[i]], ...) 中的错误：找不到对象“stim_lvis3”

原文

I have extracted time series from a few regions of interest in the brain (fMRI) and I have added pairwise correlation (Fisher-Z) values for each subject under columns corresponding the correlation between two nodes in the brain (for example: stim_lvis3, stim = stimulation site and lvis3= left visual network 3). Now, I would like to perform ANOVAs on this dataset to look at the effects and between/within group differences (3 groups x 3 timepoints). My data is already in long format.

*groups= ctbs [10 subjects x 3 timepoints], itbs = [10 subjects x 3 timepoints], and sham [10 subjects x 3 timepoints]

Any suggestions on how this can be done, given that I have 12 columns with connectivity values (stim_lvis3....stim_rpcc1). for example I have not been able to box plot the data faceted both by time and group?
How to perform a two-way mixed anova in this situation for all 12 columns for each group at a specific timepoints and then compare groups at each timepoint?

data in long format

I converted subject, time and group to factors

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

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

眉目亦如画i 2025-01-23 23:05:56

欢迎来到SO！您看起来像是新手，但将来要快速获得出色的答案，请确保以可用的格式包含示例数据（即 dput(head(myData)) 的输出）。看看：制作 R 可重现的问题。

我知道有两种在方差分析内和之间完成的方法。更容易实现的版本来自 ez 包。 jmv 包提供了更复杂的编写，但您也拥有大量的控制权。我确信 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。
对于这种方法，您无法获得事后比较；效应大小是广义的 η²。

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
# 
# 欢迎来到SO！您看起来像是新手，但将来要快速获得出色的答案，请确保以可用的格式包含示例数据（即 dput(head(myData)) 的输出） 。看看：制作 R 可重现的问题。
我知道有两种在方差分析内和之间完成的方法。更容易实现的版本来自 ez 包。 jmv 包提供了更复杂的编写，但您也拥有大量的控制权。我确信 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。

对于这种方法，您无法获得事后比较；效应大小是广义的 η²。
Mauchly's Test for Sphericity`
#       Effect         W         p p<.05
# 3       time 0.9914977 0.9541232      
# 4 group:time 0.9914977 0.9541232      
# 
# 欢迎来到SO！您看起来像是新手，但将来要快速获得出色的答案，请确保以可用的格式包含示例数据（即 dput(head(myData)) 的输出） 。看看：制作 R 可重现的问题。
我知道有两种在方差分析内和之间完成的方法。更容易实现的版本来自 ez 包。 jmv 包提供了更复杂的编写，但您也拥有大量的控制权。我确信 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。

对于这种方法，您无法获得事后比较；效应大小是广义的 η²。
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，您需要将数据旋转得更宽。您的受试者内数据需要位于单独的列中。由于时间是用stim...列中的受试者内值表示的因素，因此这就是您需要进行旋转的因素。

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 相当于科目间； 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 提供的内容之外，它还提供了每个组、时间、因变量以及每个变量的混合的事后比较；统计假设的结果：用于同质性的 Levene 和用于正态性的残差 QQ 图；以及估计的边际平均值表和这些平均值的图。

我意识到您的数据中的字段比我这里的字段多。我相信我已经为你提供了足够的基础。

我建议您确定要从数据中学习什么，并从中选择算法。

如果您有任何疑问，请告诉我。

Welcome to SO! It's looks like you're new, but in the future to get great answers quickly, make sure you include sample data in a format that useable (i.e., the output from dput(head(myData))). Check it out: making R reproducible questions.

I know of two approaches to completing within and between ANOVA. The easier to implement version is from the package ez. The package jmv offers a much more complex write-up but you have an immense amount of control, as well. I'm sure there is more to ez's version, but I haven't worked with that package very much.

I created data to somewhat simulate what you are working with.

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)

To use ez, it's a pretty straightforward implementation. Although, I couldn't find an option that allowed for multiple dependent variables. Essentially, you would either need to pivot_long or you can use jmv.
For this method, you don't get the post hoc comparisons; the effect size is the generalized η².

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
# 
# Welcome to SO! It's looks like you're new, but in the future to get great answers quickly, make sure you include sample data in a format that useable (i.e., the output from dput(head(myData))). Check it out: making R reproducible questions.
I know of two approaches to completing within and between ANOVA. The easier to implement version is from the package ez. The package jmv offers a much more complex write-up but you have an immense amount of control, as well. I'm sure there is more to ez's version, but I haven't worked with that package very much.
I created data to somewhat simulate what you are working with.
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)

To use ez, it's a pretty straightforward implementation. Although, I couldn't find an option that allowed for multiple dependent variables. Essentially, you would either need to pivot_long or you can use jmv.

For this method, you don't get the post hoc comparisons; the effect size is the generalized η².
Mauchly's Test for Sphericity`
#       Effect         W         p p<.05
# 3       time 0.9914977 0.9541232      
# 4 group:time 0.9914977 0.9541232      
# 
# Welcome to SO! It's looks like you're new, but in the future to get great answers quickly, make sure you include sample data in a format that useable (i.e., the output from dput(head(myData))). Check it out: making R reproducible questions.
I know of two approaches to completing within and between ANOVA. The easier to implement version is from the package ez. The package jmv offers a much more complex write-up but you have an immense amount of control, as well. I'm sure there is more to ez's version, but I haven't worked with that package very much.
I created data to somewhat simulate what you are working with.
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)

To use ez, it's a pretty straightforward implementation. Although, I couldn't find an option that allowed for multiple dependent variables. Essentially, you would either need to pivot_long or you can use jmv.

For this method, you don't get the post hoc comparisons; the effect size is the generalized η².
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

Now to use jmv, you will need to pivot the data wider. Your within-subject data needs to be in separate columns. Since time is the factor that represents with the within-subject values in the stim... columns, that's what you need to pivot.

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

Now you're ready for jmv. bs equates to between-subjects; rm equates to repeated measures.

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

If you print fit, you're going to get a ton of information. In addition to the content provided by ezANOVA, this provides the post hoc comparisons of each group, time, dependent variable, and the intermix of each; the results of the statistical assumptions: Levene's for homogeneity and a QQ plot of the residuals for normality; and the estimated marginal means table and a plot of those means.

I realize that you have more fields in your data than what I have here. I believe I've provided enough for you to build off of.

I suggest that you determine what you want to learn from the data and choose the algorithm from there.

If you have any questions, let me know.

回复收藏 0 原文

~没有更多了~