重复测量 - 测量处理之间差异的正确代码（随着时间的推移）

Question

对于一个研究项目，我需要运行方差分析检验以查看某些处理之间差异的统计显着性。 该实验包括将一些细菌接种在不同的试管中，这些试管含有不同浓度的不同处理剂。 我的因变量是在分光光度计上测量的光密度 660 的值，我在不同时间随时间测量了 13 次 OD。

这是数据集，我给你所有的数据集，它不是那么大：

od34_stat1 <- data.frame(
          OD = c(0.032667,0.09,0.157,0.184,0.345667,
                 0.4445,0.47725,0.53925,0.74,0.750667,0.859167,0.880333,
                 0.8275,0.034667,0.0935,0.146,0.1725,0.522167,0.5865,0.71075,
                 0.69875,0.927,0.929667,1.063167,1.037333,0.973,0.031167,
                 0.1045,0.139,0.1665,0.425667,0.523,0.69875,0.80575,
                 1.0435,0.994667,1.085667,1.215333,1.1145,0.034667,0.1085,
                 0.1285,0.1645,0.349667,0.474,0.74075,0.78125,1.0815,
                 0.937167,1.045667,1.104333,0.9555,0.028167,0.065,0.13,0.1715,
                 0.331667,0.4015,0.45775,0.54425,0.811,0.739167,0.797167,
                 0.773333,0.6905,0.021167,0.0835,0.131,0.1585,0.279167,
                 0.384,0.40225,0.46975,0.646,0.625667,0.684667,0.701333,
                 0.5885,0.015667,0.0655,0.086,0.12,0.191667,0.261,0.29875,
                 0.35825,0.446,0.411167,0.364667,0.369333,0.31),
   Treatment = as.factor(c("0_CNTRL","0_CNTRL",
                           "0_CNTRL","0_CNTRL","0_CNTRL","0_CNTRL","0_CNTRL",
                           "0_CNTRL","0_CNTRL","0_CNTRL","0_CNTRL",
                           "0_CNTRL","0_CNTRL","10_TOX","10_TOX","10_TOX","10_TOX",
                           "10_TOX","10_TOX","10_TOX","10_TOX","10_TOX",
                           "10_TOX","10_TOX","10_TOX","10_TOX","25_TOX",
                           "25_TOX","25_TOX","25_TOX","25_TOX","25_TOX","25_TOX",
                           "25_TOX","25_TOX","25_TOX","25_TOX","25_TOX",
                           "25_TOX","50_TOX","50_TOX","50_TOX","50_TOX",
                           "50_TOX","50_TOX","50_TOX","50_TOX","50_TOX",
                           "50_TOX","50_TOX","50_TOX","50_TOX","10_CNTRL",
                           "10_CNTRL","10_CNTRL","10_CNTRL","10_CNTRL","10_CNTRL",
                           "10_CNTRL","10_CNTRL","10_CNTRL","10_CNTRL",
                           "10_CNTRL","10_CNTRL","10_CNTRL","25_CNTRL","25_CNTRL",
                           "25_CNTRL","25_CNTRL","25_CNTRL","25_CNTRL",
                           "25_CNTRL","25_CNTRL","25_CNTRL","25_CNTRL",
                           "25_CNTRL","25_CNTRL","25_CNTRL","50_CNTRL","50_CNTRL",
                           "50_CNTRL","50_CNTRL","50_CNTRL","50_CNTRL",
                           "50_CNTRL","50_CNTRL","50_CNTRL","50_CNTRL","50_CNTRL",
                           "50_CNTRL","50_CNTRL")),
        Time = as.factor(c("0","2","4","6",
                           "70","94","478","496","568","616","736","784",
                           "808","0","2","4","6","70","94","478","496",
                           "568","616","736","784","808","0","2","4","6",
                           "70","94","478","496","568","616","736","784",
                           "808","0","2","4","6","70","94","478","496",
                           "568","616","736","784","808","0","2","4",
                           "6","70","94","478","496","568","616","736",
                           "784","808","0","2","4","6","70","94","478",
                           "496","568","616","736","784","808","0","2","4",
                           "6","70","94","478","496","568","616","736",
                           "784","808"))
)

所以，我尝试做的是重复测量方差分析，考虑到我随时间测量了 OD，时间是我的重复测量因子（？）。

我需要查看治疗组之间是否存在统计学上的显着差异（例如 0_CNTRL 和 25_TOX 之间是否存在显着差异？）。 最初我找到了一个代码，它在重复测量中正确执行方差分析，但它向我展示了时间点之间的差异：然后它告诉我时间 4 和时间 6 等之间是否存在差异，但这不是我要问的问题需要，最重要的是结果太分散了。

这是原始代码（我遵循了本指南：https://www.datanovia.com/en/lessons/repeated-measures-anova-in-r/#one-way-repeated-measures-anova）：

library(tidyverse)
library(ggpubr)
library(rstatix)
library(ggplot2)

##Factors
od34_stat1$Treatment <- as.factor(od34_stat1$Treatment)
od34_stat1$Time <- as.factor(od34_stat1$Time)

#Interactionplot - Boxplot
bxp34 <- ggboxplot(od34_stat1, x = "Time", y = "OD", add = "point")
bxp34

##Check assumptions: Outliers
od34_stat1 %>%
  group_by(Time) %>%
  identify_outliers(OD)

##Check assumptions: Normality
od34_stat1 %>%
  group_by(Time) %>%
 shapiro_test(OD)
#OR
ggqqplot(od34_stat1, "OD", facet.by = "Time")

#Computing One-Way repeated measure ANOVA
od34.aov <- anova_test(data = od34_stat1, dv = OD, wid = Treatment, within = Time)
get_anova_table(od34.aov)

# Pairwise comparisons
od34.pwc <- od34_stat1 %>%
  pairwise_t_test(
    OD ~ Time, paired = TRUE,
    p.adjust.method = "bonferroni"
    )
od34.pwc

##Creating Report
od34.pwc <- od34.pwc %>% add_xy_position(x = "Time")
bxp34 + 
  stat_pvalue_manual(od34.pwc) +
  labs(
    subtitle = get_test_label(od34.aov, detailed = TRUE),
    caption = get_pwc_label(od34.pwc)
  )

好的。这是我的问题，现在输出是“时间”因素。然而，该指南使用的数据集只有 3 次因变量测量，而我测量了 13 次。此外，我认为该指南的目的恰恰是要查看随时间变化的差异，而我的目的是查看随时间测量 OD 测量值的处理之间的差异..

因此，作为 rstudio 菜鸟，我的想法是将代码：“Time”更改为“Treatment”。这样我的输出就是我需要的。 我担心的是，通过改变这些因素，结果很清楚，但没有逻辑意义。

审核代码：

#Interactionplot - Boxplot
bxp34_1 <- ggboxplot(od34_stat1, x = "Treatment", y = "OD", add = "point")
bxp34_1

##Check assumptions: Outliers
od34_stat1 %>%
  group_by(Time) %>%
  identify_outliers(OD)

##Check assumptions: Normality
od34_stat1 %>%
  group_by(Treatment) %>%
 shapiro_test(OD)
#OR
ggqqplot(od34_stat1, "OD", facet.by = "Treatment")

#Computing One-Way repeated measure ANOVA
od34.aov_1 <- anova_test(data = od34_stat1, dv = OD, wid = Time, within = Treatment)
get_anova_table(od34.aov_1)

# Pairwise comparisons
od34.pwc_1 <- od34_stat1 %>%
  pairwise_t_test(
    OD ~ Treatment, paired = TRUE,
    p.adjust.method = "bonferroni"
    )
od34.pwc_1

##Creating Report
od34.pwc_1 <- od34.pwc_1 %>% add_xy_position(x = "Treatment")
bxp34_1 + 
  stat_pvalue_manual(od34.pwc_1) +
  labs(
    subtitle = get_test_label(od34.aov_1, detailed = TRUE),
    caption = get_pwc_label(od34.pwc_1)
  )

这样我的图形输出 (od34.pwc_1) 允许我解释治疗之间差异的统计显着性。

我希望我已经正确地总结了所有的疑问。你怎么认为？这样做对吗？ 如果不正确，您会推荐什么来分析和可视化这些治疗之间的差异？