R中的aov()错误术语:bw错误(id)和错误(id/timevar)规范有什么区别？

Question

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)) 

为什么会这样？哪种变体是正确的？

Answer 1

这是一篇博客文章,它将帮助分解"经典方差分析中的随机效应"部分中的每个含义.

在博客中,这里是对Error术语中"划分"的概述.

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但随后timevar在timevar嵌套的id水平与之间进行修复

aov(depvar~timevar+Error(id))

这只是id随机效应,对其他变量没有约束.

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

这可能也证明是有用的,这是一些代码分析方差,它有一些关于学习ANOVA的建议.

Answer 2

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

\n\n

请注意，将变量作为随机效应包含在内的方法不止一种。您还可以使用aov(depvar~timevar+Error(id*timevar))或aov(depvar~timevar+Error(id + timevar))。其中每一个都意味着完全不同的东西，但它可能会令人困惑，因为由于数据本身的限制，当应用于相同的数据集时，它们通常会给出相似的结果。

\n\n

/中使用的斜杠aov()表示嵌套。当您使用 时/，R 会自动将其扩展为底部变量的主效应加上底部和顶部之间的相互作用。例如，A/B自动扩展为A + A:B. 这类似于A*B自动扩展为 的方式A + B + A:B，但是通过嵌套，嵌套中的变量永远不会出现在其嵌套之外（即，不会产生主效应）B方式，但通过嵌套，嵌套中的变量永远不会出现在其嵌套之外（即，本身

\n\n

您可以在输出中看到这种扩展：

\n\n

> summary(aov(var~time1+Error(id / time1)))\n\nError: id\n          Df Sum Sq Mean Sq F value Pr(>F)\nResiduals  1  52.24   52.24               \n\nError: id:time1\n      Df Sum Sq Mean Sq\ntime1  1   4291    4291\n\nError: Within\n           Df Sum Sq Mean Sq F value  Pr(>F)   \ntime1       1   1239  1238.7   10.19 0.00167 **\nResiduals 176  21399   121.6                   \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n

\n\n

这些Error术语表示随机效应。请注意，您会得到一个表示 的主要效果，id因为它是嵌套的基础，还有一个表示id和之间的交互作用time1，因为time1嵌套在其中id（您还会得到以下Error效果：Within其中是模型的基本残差项，即个体观察本身的随机效应）。

\n\n

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

\n\n

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

\n\n

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

\n\n

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

\n\n

我假设 指的id是单个参与者，并且 和 上的 9 次测量time1是time2对每个参与者的受试者内测试（即每个参与者在 上测量了 9 次var，因此这是重复测量设计）。

\n\n

具体来说，假设var是某个问题解决任务的分数， 和time1分别time2是允许参与者研究问题的分钟数和完成问题的时间量。自从time1和time2是交叉的，每个参与者在每种情况的组合下完成任务 9 次。

\n\n

> table(time1, time2)\n     time2\ntime1  8 15 20\n    1 20 20 20\n    2 20 20 20\n    3 20 20 20\n> table(time1, time2, id)\n, , id = 1\n\n     time2\ntime1 8 15 20\n    1 1  1  1\n    2 1  1  1\n    3 1  1  1\n\n, , id = 2\n\n     time2\ntime1 8 15 20\n    1 1  1  1\n    2 1  1  1\n    3 1  1  1\n(output truncated)\n

\n\n

参与者进行分组测试，一半参与者在第一组，另一半在第二组。也许研究是在教室进行的，第一组是一个班级，第二组是第二个班级。也许，群体身份实际上并不是一个感兴趣的变量，但我们不应该将其排除在模型之外，因为群体之间的差异可能会导致一些令人讨厌的方差。例如，也许第一间教室有更好的照明，使第一组的所有成员比第二组的成员更有机会在拼图上取得好成绩。

\n\n

分数、ID 和组都应该是随机效应，时间 1 和时间 2 应该是固定效应（请注意，如果您对模型有不同的想法，则对于相同的数据，这可能会有所不同；例如，您可能希望根据您的情况将组视为固定效应）研究问题）。

\n\n

给定该模型，这将是该模型最完整的版本，使用aov()：

\n\n

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

\n\n

这是输出：

\n\n

> summary(aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df))\n\nError: group\n          Df Sum Sq Mean Sq F value Pr(>F)\nResiduals  1  771.7   771.7               \n\nError: group:id\n          Df Sum Sq Mean Sq F value Pr(>F)\nResiduals 18  243.8   13.55               \n\nError: group:id:time1\n          Df Sum Sq Mean Sq F value Pr(>F)    \ntime1      2   7141    3571   181.6 <2e-16 ***\nResiduals 38    747      20                   \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n\nError: group:id:time2\n          Df Sum Sq Mean Sq F value Pr(>F)    \ntime2      2  16353    8176   434.6 <2e-16 ***\nResiduals 38    715      19                   \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n\nError: group:id:time1:time2\n            Df Sum Sq Mean Sq F value  Pr(>F)   \ntime1:time2  4  214.5   53.63   5.131 0.00103 **\nResiduals   76  794.3   10.45                   \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\nWarning message:\nIn aov(var ~ time1 * time2 + Error(group/id/(time1 * time2)), data = df) :\n  Error() model is singular\n

\n\n

（除了上面的链接之外，这里还有一些关于随机效应与固定效应的额外指导）

\n