logo标签列表

R中的aov()函数错误项：Error(id)和Error(id/timevar)规范有什么区别？

ranova
17
aov(depvar~timevar+Error(id)) 和 aov(depvar~timevar+Error(id/timevar)) 公式规范有何区别？这两种变体会产生略微不同的结果。
同样的问题曾在此处提问过：https://stats.stackexchange.com/questions/60108/how-to-write-the-error-term-in-repeated-measures-anova-in-r不过，我想用一个更恰当的例子再次重复一遍。
下面是我创建的一个例子：
var=rep(NA,180)
id=rep(1:20,each=180/20)
group=rep(rep(1:2,each=9),180/(9*2))
time1=rep(rep(1:3,each=3),180/(3*3))
time2=rep(c(8,15,20),180/3)

var[group==1&time1==1&time2==8]=runif(10,105,115)
var[group==2&time1==1&time2==8]=runif(10,105,115)
var[group==1&time1==1&time2==15]=runif(10,95,105)
var[group==2&time1==1&time2==15]=runif(10,95,105)
var[group==1&time1==1&time2==20]=runif(10,85,95)
var[group==2&time1==1&time2==20]=runif(10,85,95)

var[group==1&time1==2&time2==8]=runif(10,95,105)
var[group==2&time1==2&time2==8]=runif(10,95,105)
var[group==1&time1==2&time2==15]=runif(10,85,95)
var[group==2&time1==2&time2==15]=runif(10,75,85)
var[group==1&time1==2&time2==20]=runif(10,75,85)
var[group==2&time1==2&time2==20]=runif(10,65,75)

var[group==1&time1==3&time2==8]=runif(10,95,105)
var[group==2&time1==3&time2==8]=runif(10,95,105)
var[group==1&time1==3&time2==15]=runif(10,85,95)
var[group==2&time1==3&time2==15]=runif(10,75,85)
var[group==1&time1==3&time2==20]=runif(10,75,85)
var[group==2&time1==3&time2==20]=runif(10,65,75)

df=data.frame(id,var,group,time1,time2)
df$id=factor(df$id)
df$group=factor(df$group)
df$time1=factor(df$time1)
df$time2=factor(df$time2)

在此执行aov()会根据Error()项的不同规定获得稍微不同的结果：

仅针对一个时间项：

> summary(aov(var~time1+Error(id),data=df))
Error: id
      Df Sum Sq Mean Sq F value Pr(>F)
      Residuals 19  958.4   50.44               
Error: Within
       Df Sum Sq Mean Sq F value   Pr(>F)    
       time1       2   7538    3769   30.41 6.72e-12 ***
       Residuals 158  19584     124         

> summary(aov(var~time1+Error(id/time1),data=df))
Error: id
      Df Sum Sq Mean Sq F value Pr(>F)
      Residuals 19  958.4   50.44               
Error: id:time1
      Df Sum Sq Mean Sq F value Pr(>F)    
      time1      2   7538    3769   211.5 <2e-16 ***
      Residuals 38    677      18                   
      ---
     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
       Df Sum Sq Mean Sq F value Pr(>F)
       Residuals 120  18907   157.6

或者对于两个时间术语（为了节省空间，请不要在此处输入输出，您可以自行检查）：

summary(aov(var~group*time1*time2+Error(id/(group*time1*time2)),data=df)) 
summary(aov(var~group*time1*time2+Error(id),data=df))

为什么会发生这种情况？哪个变量是正确的？
2
这个问题真正属于stats.stackexchange.com。 - Chris
2个回答
14

以下是一篇博客文章，将有助于在“经典ANOVA中的随机效应”部分中解释每个术语的含义。

从该博客中，以下是关于“误差”项中“dividing”的摘要。

aov(Y ~ Error(A), data=d)               # Lone random effect
aov(Y ~ B + Error(A/B), data=d)         # A random, B fixed, B nested within A
aov(Y ~ (B*X) + Error(A/(B*X)), data=d) # B and X interact within levels of A

所以根据你的问题，

aov(depvar~timevar+Error(id/timevar))

这意味着您对 id 产生了随机效应，但是使用在 id 层次内嵌套的 timevar 来固定 timevar。

aov(depvar~timevar+Error(id))

这只是将id作为随机效应，没有对其他变量进行限制。

来源：http://conjugateprior.org/2013/01/formulae-in-r-anova/

这个链接也许会有用，其中有一些代码介绍方差分析，并提供了学习ANOVA的一些建议。

6

aov(depvar~timevar+Error(id))和aov(depvar~timevar+Error(id/timevar))之间的区别在于是否将timevar作为随机效应包含在内。

请注意，有多种方法可以将变量作为随机效应包含在内。你还可以使用aov(depvar~timevar+Error(id*timevar))或者aov(depvar~timevar+Error(id + timevar))。虽然这些方式各有不同的含义，但是由于数据本身的限制，它们通常在应用于相同的数据集时会给出类似的结果，这可能会让人感到困惑。

在 aov() 中使用的斜杠 / 表示 嵌套。当您使用 / 时，R 会自动将其扩展为底部变量的主效应和底部与顶部之间的交互作用。例如，A/B 自动扩展为 A + A:B。这类似于 A*B 如何自动扩展为 A + B + A:B，但是在嵌套中，嵌套中的变量永远不会出现在其嵌套之外（即不能有 B 的主效应）。您可以在输出中看到此扩展发生的情况：
> summary(aov(var~time1+Error(id / time1)))

Error: id
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  1  52.24   52.24               

Error: id:time1
      Df Sum Sq Mean Sq
time1  1   4291    4291

Error: Within
           Df Sum Sq Mean Sq F value  Pr(>F)   
time1       1   1239  1238.7   10.19 0.00167 **
Residuals 176  21399   121.6                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error术语表示随机效应。请注意，对于id的主要效应，您会得到一个，因为它是嵌套的基础，对于id和time1之间的交互作用，您也会得到一个，因为time1嵌套在id中（您还会获得Within的Error效应，这是模型的基本残差项，即个体观察的随机效应）。

那么，您的数据的正确方法是什么？

这取决于1）您的数据实际上的结构以及2）您打算运行的模型。注意：没有明确的测试可以在数据上运行以确定结构或正确的模型；这是一种思考练习而不是计算练习。

在您提供的示例模型中，您有一个结果var，然后是看起来像分组变量group和id，以及两个时间变量time1和time2。每个id仅在1个组中，而不是跨越两个组，这表明id嵌套在组内。

> table(group, id)
     id
group 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
    1 9 0 9 0 9 0 9 0 9  0  9  0  9  0  9  0  9  0  9  0
    2 0 9 0 9 0 9 0 9 0  9  0  9  0  9  0  9  0  9  0  9

我假设id指的是单个参与者，而time1和time2上的9项测量是在每个参与者内进行的测试（即每个参与者在var上被测量了9次，因此这是一种重复测量设计）。
具体来说，假设var是一项问题解决任务的得分，而time1和time2分别是参与者被允许研究问题的时间和完成问题所需的时间。由于time1和time2是交叉的，每个参与者在各种情况下都要完成任务9次。
> table(time1, time2)
     time2
time1  8 15 20
    1 20 20 20
    2 20 20 20
    3 20 20 20
> table(time1, time2, id)
, , id = 1

     time2
time1 8 15 20
    1 1  1  1
    2 1  1  1
    3 1  1  1

, , id = 2

     time2
time1 8 15 20
    1 1  1  1
    2 1  1  1
    3 1  1  1
(output truncated)

参与者分为两组进行测试，一半参与者在第一组，另一半在第二组。也许这项研究是在教室里进行的，第一组是一个班级，第二组是第二个班级。可能实际上并不关心组别身份，但我们不应该将其从模型中排除，因为组别之间的差异可能会导致某些麻烦的方差。例如，也许第一个教室有更好的照明，给第一组成员比第二组成员更好的完成难题的机会。
分数、ID和组别都应该是随机效应，而time1和time2应该是固定效应（请注意，如果对同一数据有不同的模型想法，则此项可能会有所变化；例如，根据你的研究问题，你可能希望将组别视为固定效应）。
基于该模型，使用aov()，这将是最完整的模型版本。
aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df)

这是输出结果：
> summary(aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df))

Error: group
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  1  771.7   771.7               

Error: group:id
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 18  243.8   13.55               

Error: group:id:time1
          Df Sum Sq Mean Sq F value Pr(>F)    
time1      2   7141    3571   181.6 <2e-16 ***
Residuals 38    747      20                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: group:id:time2
          Df Sum Sq Mean Sq F value Pr(>F)    
time2      2  16353    8176   434.6 <2e-16 ***
Residuals 38    715      19                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: group:id:time1:time2
            Df Sum Sq Mean Sq F value  Pr(>F)   
time1:time2  4  214.5   53.63   5.131 0.00103 **
Residuals   76  794.3   10.45                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In aov(var ~ time1 * time2 + Error(group/id/(time1 * time2)), data = df) :
  Error() model is singular

除了以上链接，这里还有一些关于随机效应 vs. 固定效应的额外指导。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接