您正确地注意到，通常的 ANOVA 无法处理这种非常异方差且高度依赖的时间点数据。但是，其他重复测量或多变量过程也会失败，因为您只有 6 个重复但有 10 个时间点，即所谓的高维数据。最后证明，不能对这样的数据进行精确的测试。

但是，有一个很好的近似测试，它是 Huynh-Feldt 程序，请参阅此技术报告以获取参考资料和对您可能使用的可能不等协方差矩阵的概括。它甚至是为您非常小的数据集定义的。没有关于方差和协方差的假设。它适用于正态分布的数据，但我认为你可以接受。

您可以测试曲线是否不同 (set T = diag(10)) 以及治疗和时间点之间是否存在交互 (set T = diag(10) - 1/10)。对于额外的R-package 来说太简单了。检查这个：

Tmat <- diag(10) # or Tmat <- diag(10)-1/10
Y1 <- t(matrix(my.data[my.data$treatment=="B","value"],10,3)) %*% Tmat
    Y2 <- t(matrix(my.data[my.data$treatment=="A","value"],10,3)) %*% Tmat

Mean1 <- apply(Y1,2,mean)
Mean2 <- apply(Y2,2,mean)
S1 <- cov(Y1)
S2 <- cov(Y2)

F <- sum((Mean1-Mean2)**2) / (sum(diag(S1/3 + S2/3))) # see eq. 3.18
trS1sq <- 3/2 * ( (sum(diag(S1)))**2 - 2/3 * sum(S1**2)) # see eq. 3.26
trS2sq <- 3/2 * ( (sum(diag(S2)))**2 - 2/3 * sum(S2**2)) 
trS1S1 <- 3/2 * ( sum(S1**2) - (sum(diag(S1)))**2/2)  # 3.27
trS2S2 <- 3/2 * ( sum(S2**2) - (sum(diag(S2)))**2/2)  # 3.27

f <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1+trS2S2+2*sum(S1*S2))  # see p. 16
f0 <- (trS1sq + trS2sq + 2*sum(diag(S1))*sum(diag(S2))) / 
    (trS1S1/2 + trS2S2/2)

p.value <- 1-pf(F,f,f0)

我相信您的问题是纵向数据分析的教科书案例，为此开发了广泛的统计方法。特别是，分析响应配置文件是我在这种情况下会采用的一种方法。

为了激发这一点，我将提供一些来自GM Fitzmaurice的应用纵向分析的背景：

1.2 纵向和聚类数据

纵向研究的决定性特征是对同一个人的测量随着时间的推移而重复进行，从而可以直接研究随时间的变化。纵向研究的主要目标是描述响应随时间的变化以及影响变化的因素。(...)

纵向数据的一个显着特征是它们是 聚集的。在纵向研究中，集群由在不同场合从单个个体获得的重复测量值组成。集群内的观察通常会表现出正相关，并且必须在分析中考虑这种相关性。纵向数据也有时间顺序；集群内的第一次测量必须在第二次测量之前，依此类推。重复测量的顺序对分析具有重要意义。

以及关于分析响应配置文件的选择：

当存在单个分类协变量（可能表示不同的治疗或暴露组）并且无法指定组间反应谱差异的特定先验模式时，分析反应谱的方法很有吸引力。当在相同的场合序列中获得重复测量时，可以通过每个场合的估计平均响应来总结数据，并按组因子的水平进行分层。在组因子的任何给定水平上，随时间变化的均值序列称为均值响应曲线。

该nlme R库包含一个glm评估此类模型的函数，为您提供 Anova 的纵向“等效”：

library(nlme)

# Model assuming the same variance for each time point
gls.fit <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit)

# Model allowing for different variance structure for each time point 
gls.fit.diff.var <- 
  gls(value ~ factor(time) + treatment, 
      data = my.data,
      corr = corSymm(form = ~ 1 | object),
      weights = varIdent(form = ~ 1 | factor(time)),
      control = glsControl(tolerance = 0.01, msTol = 0.01, 
                           maxIter = 1000000, msMaxIter = 1000000))
summary(gls.fit.diff.var)

不幸的是，在这种情况下，模型估计没有覆盖（即使我control为了更方便更改了参数）。恐怕您的数据集中的案例太少，无法估计参数。