如何在 R 中为嵌套在主题设计中获得球形度？

机器算法验证 r 方差分析 球形度
2022-03-15 11:17:28

感谢您的阅读。我试图获得纯粹在主题设计中的球形值。我一直无法使用 ezANOVA 或 Anova()。如果我添加一个主体之间的因素，Anova 就可以工作，但我无法获得纯粹在主体设计中的球形度。有什么建议吗？

我已经阅读了 R 时事通讯、狐狸章节附录、EZanova 以及我可以在网上找到的任何内容。

我原来的方差分析

anova(aov(resp ~ sucrose*citral, random =~1 | subject, data = p12bl, subset = exps==1)) 
anova(aov(resp ~ sucrose*citral, random =~1 | subject/sucrose*citral, data = p12bl, subset = exps==1))

> str(subset(p12bl, exps==1))
'data.frame':   192 obs. of  12 variables:
 $ exps     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Order    : int  1 1 1 1 1 1 1 1 1 1 ...
 $ threshold: Factor w/ 2 levels " Suprathreshold",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ SET      : Factor w/ 2 levels " A"," B": 1 1 1 1 1 1 1 1 1 1 ...
 $ subject  : Factor w/ 16 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ stim     : chr  "S1C1" "S1C1" "S1C1" "S1C1" ...
 $ resp     : num  6.01 5.63 0 2.57 6.81 ...
 $ id       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ X1       : Factor w/ 1 level "S": 1 1 1 1 1 1 1 1 1 1 ...
 $ sucrose  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ X3       : Factor w/ 1 level "C": 1 1 1 1 1 1 1 1 1 1 ...
 $ citral   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

subset(p12b,exps==1)
   exps Order       threshold SET observ S1C1 S1C2 S1C3 S1C4 S2C1 S2C2 S2C3 S2C4 S3C1 S3C2 S3C3 S3C4 S4C1 S4C2 S4C3 S4C4
1     1     1  Suprathreshold   A      1  6.0  7.1  7.5  8.6 15.0 15.4 15.0 13.1 16.9   13 13.1 16.5   24   16   21   20
2     1     1  Suprathreshold   A      2  5.6  0.8  4.0  5.6  5.6 11.3 12.9 14.5 18.5   15 12.9 14.5   24   26   29   28
3     1     1  Suprathreshold   A      3  0.0  0.0  1.7  0.0  5.0  8.4  8.4  5.0 11.7   20 18.5 16.8   29   37   37   30
4     1     1  Suprathreshold   A      4  2.6  3.3  9.1 16.3  5.4 10.0  9.6 16.8 13.5   12 22.2 23.1   19   20   22   23
5     1     1  Suprathreshold   A      5  6.8  5.3 15.4 14.5 11.5  8.3 14.5 14.2  8.9   17 11.2 15.1   24   23   19   19
6     1     1  Suprathreshold   A      6  2.6  2.8  2.6  5.2 13.4 15.6 13.7 13.0 13.7   15 16.0 18.9   22   24   25   25
7     1     1  Suprathreshold   A      7  1.3  5.8 10.2  9.8 11.9 12.3 17.7 16.7 11.4   19 19.2 21.1   16   19   18   19
8     1     1  Suprathreshold   A      8  2.0  5.6  3.9  2.0  4.9  5.2  7.5  4.9 20.2   21  8.2  9.5   30   26   32   45
9     1     1  Suprathreshold   A      9  9.4 11.3 11.7 12.1 14.7 13.8 12.6 14.9 15.2   15 15.9 13.9   17   18   15   18
10    1     1  Suprathreshold   A     10  4.5 17.8 18.5 21.6  5.8 10.9 17.0 20.2  6.6   10 17.8 18.7   12   12   16   19
11    1     1  Suprathreshold   A     11  9.8 13.0 16.1 18.0 10.5 11.6 15.4 17.3 10.1   14 15.2 16.7   13   15   15   17
12    1     1  Suprathreshold   A     12  9.6 10.4 13.3 11.3 12.1 12.6 13.6 13.6 14.9   16 15.1 16.3   16   18   18   17

样本输出

ezANOVA( data = subset(p12bl, exps==1)  , dv= .(resp), sid = .(observ), within = .(sucrose,citral), between = NULL, collapse_within = FALSE)
Note: model has only an intercept; equivalent type-III tests substituted.
$ANOVA
          Effect DFn DFd  SSn  SSd     F       p p<.05   pes
1        sucrose   3  33 4953 3263 16.70 9.0e-07     * 0.603
2         citral   3  33  410  553  8.16 3.3e-04     * 0.426
3 sucrose:citral   9  99   56  791  0.77 6.4e-01       0.066

警告消息： 1：您已从分析中删除了一个或多个 S。重构方差分析的“观察”。2：Anova() 的 Ss 太少，恢复到 aov()。请参阅 ezANOVA 帮助的“警告”部分。

4个回答

尝试：

library(ez)
ezANOVA(data=subset(p12bl, exps==1),
  within=.(sucrose, citral),
  wid=.(subject),
  dv=.(resp)
  )

和等效的 aov 命令，减去球形等，是：

aov(resp ~ sucrose*citral + Error(subject/(sucrose*citral)), 
    data= subset(p12bl, exps==1))

这是直接从汽车使用 Anova 的等价物：

library(car)
df1<-read.table("clipboard", header=T) #From copying data in the question above
sucrose<-factor(rep(c(1:4), each=4))
citral<-factor(rep(c(1:4), 4))
idata<-data.frame(sucrose,citral)

mod<-lm(cbind(S1C1, S1C2, S1C3, S1C4, S2C1, S2C2, S2C3, S2C4, 
        S3C1, S3C2, S3C3, S3C4, S4C1, S4C2, S4C3, S4C4)~1, data=df1)
av.mod<-Anova(mod, idata=idata, idesign=~sucrose*citral)
summary(av.mod)

你试过car约翰福克斯的包裹吗？它包括在Anova()处理实验设计时非常有用的功能。它应该在 Greenhouse-Geisser 校正和 Huynh-Feldt 校正之后为您提供校正的 p 值。如果您想知道如何使用它，我可以发布一个快速的 R 示例。

此外，还有一个很好的教程，介绍了使用 R 进行重复测量和混合效应模型进行心理学实验和问卷调查；见第 6.10 节关于球形度。

作为旁注，Mauchly 的球形检验在 中可用，但如果我没记错的话mauchly.test()，它不适用于对象。aov2007 年 10 月的R 时事通讯包括对该主题的简要描述。

我通常建议使用现代混合建模方法完全避免这些类型的球形度测试。如果您不使用少数主题，这将为您在建模适当的协方差结构时提供很大的灵活性，从而在必要时将您从严格的球形假设中解脱出来。我从输出中推断出str您有 16 个主题，每个主题有 12 个观察值（我假设平衡 b/c 您正在使用经典矩量法工具），这应该是足够的数据来拟合具有结构化协方差矩阵的混合模型（受限）最大值可能性。

在不接近您的数据的情况下，我无法提供具体的模型建议，但是在 R 中开始的一个地方是用(after ) 替换aov您的模型规范。这会起作用的原因是您错误地提供了一个-style random 参数（正如@Matt Albrecht 指出的那样，一个术语本来是合适的）。在 nlme 中，将随机参数设置为 并且没有或参数，您正在为每个组指定一个随机截距，这意味着组内的响应协方差是lmelibrary(nlme)nlmeaovError~ 1|<your grouping structure>correlationweightZGZ' + R = 1G1'+ \sigma^2 I==> 复合对称性，组间方差非对角线和组间方差 + 组内方差对角线 ==> 球形结构。从那里您可以开始探索（例如使用内置的图形方法）、建模和测试（例如比较信息标准或使用 LRT 用于嵌套模型）各种形式的非球形。建模组件的一些工具是：

  • 使用weights参数对组内或组间的非恒定方差（对角线）建模（例如，蔗糖水平之间的误差方差变化）。
  • 使用correlation参数对组内的非常量协方差（非对角线）建模（例如，组内的结构，其中残差在时间上更接近（例如 AR1 结构）或空间（例如球形结构）更相似）。
  • 通过向公式|中的 LHS 添加项来建模随机斜率。random

虽然这个过程可能很复杂，有很多潜在的陷阱，但我相信它会引导你更多地思考数据生成机制，并结合仔细的图形检查（我推荐latticeb/cnlme具有出色的基于点阵的绘图方法——但ggplot有效也很好）你可能不仅对这个过程有更好的科学理解，而且还有更少的偏见和更有效的估算器来得出推论。

ez 现已更新至 2.0 版。在其他改进中，导致它无法在此示例中工作的错误已得到修复。

其它你可能感兴趣的问题
上一篇如何对转换后仍然不正常的数据进行方差分析？ 下一篇OLS 回归结果：p 值 > 0.10，如何进行？