R中的重复测量/受试者内组ANOVA

Question

我正试图使用R运行重复测量的ANOVA分析。我查看了各种网站上的各种示例，但它们似乎从来没有谈到我遇到的错误。我认为我可能误解了一些重要的东西。

我要运行的ANOVA是基于使用人类参与者的实验数据。它有一个因变量和三个自变量。所有自变量的所有水平都在所有参与者上运行，使其成为三因素重复测量/受试者内ANOVA。

我在R中运行的代码如下：

aov.output = aov(DV~ IV1 * IV2 * IV3 + Error(PARTICIPANT_ID / (IV1 * IV2 * IV3)),
                 data=fulldata)

当我运行这个程序时，会出现以下警告：

Error() model is singular

有没有什么想法，我可能做错了什么？

- vize

---

这个错误的快速搜索（这通常是一个好策略）将我带到了这个页面：http://tolstoy.newcastle.edu.au/R/help/04/10/5215.html 相关部分在这里：我认为这意味着正确的错误模型是Error(Subject/T.norm.Class)：我的猜测是WasSick是一个主题级别的观察，因此每个主题只有一个级别。当然，这就是拟合的模型。- Brian Ripley教授。/结束引用。我怀疑您已经指定了一个不正确的错误分布，但没有更多的信息，很难确定。 - richiemorrisroe

看起来你的随机效应部分过于复杂。奇异模型通常表示你试图在数据/观测不足的情况下拟合了一个过于复杂的模型。 - Gavin Simpson

顺便说一下，这个问题不适合在这个网站上提问 - 你最好去http://stats.stackexchange.com 上提问。 - Gavin Simpson

@richiemorrisroe 嗯，我已经谷歌搜索过了，但是错过了你指向的链接。对于所有参与者，自变量在所有水平上都有给定值，因此并不是每个受试者只有其中任何一个水平。因此，IV1有2个水平，IV2有5个水平，IV3有2个水平。 - vize

@Gavin Simpson，“太复杂的模型”是什么意思？另外，感谢您在stats.stackexchange上给出的提示-以后我会在那里发布类似的问题。我已经在那里看过了，但那里似乎更多关于纯统计和理论的问题，而不是关于R和脚本/编程的问题，这些问题似乎更适合stackoverflow。 - vize

Answer 1

尝试使用lme4软件包中的lmer函数，而aov函数可能不适用于此处。查看Douglas Bates的参考文献，例如http://lme4.r-forge.r-project.org/book/Ch4.pdf（其他章节也很棒，但那是重复测量章节，这是介绍：http://lme4.r-forge.r-project.org/book/Ch1.pdf）。R代码与相同位置，对于纵向数据，现在通常认为与其仅拟合OLS而不是像lme4包或nlme中的方差分量模型等更好。lme4似乎最近在受欢迎程度方面大幅领先于nlme。您可以注意到评论部分中Brian Ripley的引用帖子也建议转换为lme。

顺便说一下，一个巨大的优势是您将能够通过典型语法将每个效应的水平调整为总平均值的估计：

lmer(DV ~ 1  +IV1*IV2*IV3 +(IV1*IV2*IV3|Subject), dataset))

请注意，您的随机效应将是向量值。

- Patrick McCann

---

太好了，谢谢！现在看起来它正在工作。非常感谢你的提示，lmers似乎是正确的选择 :) - vize

Answer 2

我知道这篇文章的答案已经被选择了。但我仍然希望指出如何在拟合多次重复测量数据的aov或lmer模型时，如何指定正确的误差项/随机效应。我假设两个自变量（IVs）是固定的，并且交叉和与受试者一起，这意味着所有受试者都暴露于IVs的所有组合中。我将使用Kirk的《Experimental Design: Procedures for the Behavioral Sciences》（2013年）中的数据。

library(lme4)
library(foreign)
library(lmerTest)
library(dplyr)

file_name <- "http://www.ats.ucla.edu/stat/stata/examples/kirk/rbf33.dta" #1
d <- read.dta(file_name) %>%                                              #2
  mutate(a_f = factor(a), b_f = factor(b), s_f = factor(s))               #3

head(d)
    ##   a b s  y a_f b_f s_f
    ## 1 1 1 1 37   1   1   1
    ## 2 1 2 1 43   1   2   1
    ## 3 1 3 1 48   1   3   1
    ## 4 2 1 1 39   2   1   1
    ## 5 2 2 1 35   2   2   1

本研究对5个受试者（s）进行了2种不同的处理方式，包括节奏类型（a）和训练时长（b），每种处理方式各有3个水平。结果变量是对少数族裔的态度。在第3步中，我将a，b和s转化为因子变量a_f，b_f和s_f。设p和q分别为a_f和b_f的水平数量（均为3），n为样本量（5）。

在这个例子中，对于a_f，b_f以及它们的交互作用的检验应该分别使用自由度（dfs）为p-1=2、q-1=2和（p-1）*（q-1）=4。s_f误差项的df为(n-1)=4，而Within (s_f:a_f:b_f)误差项的df为(n-1) * (pq-1)=32。因此正确的模型应该给出这些dfs。

使用aov

现在让我们尝试使用不同的模型规范来使用aov：

aov(y ~ a_f*b_f + Error(s_f), data=d) %>% summary()         # m1

aov(y ~ a_f*b_f + Error(s_f/a_f:b_f), data=d) %>% summary() # m2

aov(y ~ a_f*b_f + Error(s_f/a_f*b_f), data=d) %>% summary() # m3

在m1中将错误简单地指定为Error(s_f)，可以得到与书中数值相匹配的正确dfs和F-比值。m2也给出了与m1相同的值，但还有臭名昭著的“警告：Error()模型是奇异的”。m3正在做一些奇怪的事情。它在m1（634.9）中进一步将Within残差分成三个错误项的残差：s_f:a_f (174.2)、s_f:b_f (173.6)和s_f:a_f:b_f (287.1)。这是错误的，因为当我们运行2因素组间方差分析时，我们不会得到三个错误项！多个错误项也与使用块级因子设计的目的不符，该设计允许我们对A、B和AB的测试使用相同的错误项，而分割区组设计则需要2个错误项。

使用lmer

如果数据是平衡的，那么lmerTest中使用的Kenward-Roger近似方法将给出精确的dfs和F-比值。

lmer(y ~ a_f*b_f + (1|s_f), data=d) %>% anova()         # mem1

lmer(y ~ a_f*b_f + (1|s_f/a_f:b_f), data=d) %>% anova() # mem2

lmer(y ~ a_f*b_f + (1|s_f/a_f*b_f), data=d) %>% anova() # mem3

lmer(y ~ a_f*b_f + (1|s_f:a_f:b_f), data=d) %>% anova() # mem4

lmer(y ~ a_f*b_f + (a_f*b_f|s_f), data=d) %>% anova()   # mem5

仅将随机效应指定为(1|s_f)即可获得正确的dfs和F-比率(mem1)。mem2-5甚至没有给出结果，可能需要估计的随机效应数量大于样本大小。

- Masato Nakazawa