0

注意：这个问题最初发布在Cross Validated上，建议应该在 StackOverflow 中提问。

我正在尝试对 3 向重复测量实验进行建模，FixedFactorA * FixedFactorB * Time[days]。没有遗漏的观察结果，但我的组 (FactorA * FactorB) 不相等（接近，但不完全平衡）。从在线阅读来看，对观察顺序很重要（由于响应均值和方差随时间变化）和不平等组的重复测量实验建模的最佳方法是使用混合模型并指定适当的协方差结构. 但是，我对混合模型的想法很陌生，我对是否使用正确的语法来建模我想要建模的东西感到困惑。

我想做一个完整的因子分析，这样我就可以检测到重要的时间*因子相互作用。例如，对于 FactorA = 1 的受试者，他们随时间的反应可能具有与 FactorA = 2 的受试者不同的斜率和/或截距。我还希望能够检查 FactorA 和 FactorB 的某些组合是否随着时间的推移具有显着不同的响应（因此是完整的三向交互项）。

从网上看，AR1似乎是一个合理的纵向数据协方差结构，所以我决定尝试一下。另外，我看到如果一个人计划比较两个不同的模型，应该使用 ML，所以我选择了这种方法，因为预计需要对模型进行微调。我的理解也是，目标是在模型选择期间最小化 AIC。

这是我在 SPSS（用于长格式数据）中尝试的日志中的代码，它产生的 AIC 为 2471：

MIXED RESPONSE BY FactorA FactorB Day
  /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1) SINGULAR(0.000000000001) HCONVERGE(0,
ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=FactorA FactorB Day FactorA*FactorB FactorA*Day FactorB*Day FactorA*FactorB*Day | SSTYPE(3)
  /METHOD=ML
  /PRINT=SOLUTION TESTCOV
  /REPEATED=Day | SUBJECT(Subject_ID) COVTYPE(AR1)

这是我在 R 中尝试的，它产生了 2156 的 AIC：

    require(nlme)

    #output error fix: https://stats.stackexchange.com/questions/40647/lme-error-iteration-limit-reached 

    ctrl <- lmeControl(opt='optim')  #I used this b/c otherwise I get the iteration limit reached error

    fit1 <- lme(RESPONSE ~ Day*FactorA*FactorB, random = ~ Day|Subject_ID, control=ctrl,
        correlation=corAR1(form=~Day), data, method="ML")

    summary(fit1)

这些是我的问题：

  1. 上面的 SPSS 代码产生了一个 AIC = 2471 的模型，而 R 代码产生了一个 AIC = 2156 的模型。是什么使模型不同的代码？

  2. 根据我上面的描述，这些模型中的任何一个都适合我要测试的内容吗？如果没有，什么是更好的方法，我将如何在两个程序中做到这一点以获得相同的结果？

编辑

另一件需要注意的事情是我没有对我的因素进行虚拟编码。我不知道这是否是任何一个软件的问题，或者 SPSS 与 R 中的内置编码是否不同。我也不知道这是否会成为我的三向交互术语的问题。

另外，当我说“因素”时，我指的是一个不变的群体或特征（如“性”）。

4

1 回答 1

0

从一个无条件模型开始，一个在级别 1 具有恒等方差-协方差结构，一个在级别 1 具有 AR(1) var-covar 结构：

unconditional.identity<-lme(RESPONSE~Day, random=~Day|Subject_ID, data=data, method='ML')
unconditional.ar1<-lme(RESPONSE~Day, random=~Day|Subject_ID, correlation=corAR1(form=~Day), data=data, method='ML')

求这个无条件模型的类内相关系数，即二级误差除以一级和二级误差之和。这在电子表格程序中可能更容易，但在 R 中：

intervals(unconditional.identity)$reStruct$Subject_ID[2]^2/(intervals(unconditional.identity)$reStruct$Subject_ID[2]^2+intervals(unconditional.identity)$sigma[2]^2)

intervals(unconditional.ar1)$reStruct$Subject_ID[2]^2/(intervals(unconditional.ar1)$reStruct$Subject_ID[2]^2+intervals(unconditional.ar1)$sigma[2]^2)

这取决于您的领域，但在教育研究中，低于 0.2（绝对低于 0.1）的 ICC 被认为不适合分层线性模型。也就是说，多元回归会更好，因为独立性的假设得到了证实。如果您的 ICC 低于您的领域的临界值，则不要使用分层纵向模型。

如果您的 ICC 对分层线性模型是可接受的，则添加具有恒等式和 AR(1) var-covar 矩阵的控制分组变量：

conditional1.identity<-lme(RESPONSE~Day+Group, random=~Day+Group|Subject_ID, data=data, method='ML')
conditional1.ar1<-lme(RESPONSE~Day+Group, random=~Day+Group|Subject_ID, correlation=corAR1(form=~Day), data=data, method='ML')

如果您的因素是时间不变的（您在Cross Validated上说过），那么您的模型会变得更大，因为时间和组嵌套在这些固定效应中：

conditional2.identity<-lme(RESPONSE~Day+Group+FactorA+FactorB+FactorA*Day+FactorB*Day+FactorA*Group+FactorB*Group+FactorB, random=~Day+Group|Subject_ID, data=data, method='ML')
conditional2.ar1<-lme(Day+Group+FactorA+FactorB+FactorA*Day+FactorB*Day+FactorA*Group+FactorB*Group+FactorB, random=~Day+Group|Subject_ID, correlation=corAR1(form=~Day), data=data, method='ML')

您可以获得系数的置信区间intervals()或 p 值summary()。请记住，lme以标准偏差格式报告错误术语。

我不知道你的研究领域，所以我不能说你的三向交互效应是否具有理论意义。但是此时您的模型变得非常密集。你估计的参数越多，模型在比较时的自由度就越大，所以统计显着性就会有偏差。如果你真的对三向交互效应感兴趣，我建议你考虑一下这种交互的理论意义，以及如果确实发生了这种交互意味着什么。尽管如此，您可以通过将其添加到上面的代码中来估计它：

conditional3.identity<-lme(RESPONSE~Day+Group+FactorA+FactorB+FactorA*Day+FactorB*Day+FactorA*Group+FactorB*Group+FactorB+Day*FactorA*FactorB, random=~Day+Group|Subject_ID, data=data, method='ML')
conditional3.ar1<-lme(Day+Group+FactorA+FactorB+FactorA*Day+FactorB*Day+FactorA*Group+FactorB*Group+FactorB+Day*FactorA*FactorB, random=~Day+Group|Subject_ID, correlation=corAR1(form=~Day), data=data, method='ML')

最后，比较嵌套模型：

anova(unconditional.identity,conditional1.identity,conditional2.identity,conditional3.identity)

anova(unconditional.ar1,conditional1.ar1,conditional2.ar1,conditional3.ar1)

就像我说的，你估计的参数越多，你的统计显着性就越有偏差：即，更多的参数=更多的自由度=更少的统计显着模型的机会。

但是，关于多级模型的最佳部分是比较效果大小，因此您根本不必担心 p 值。效应量采用“解释的方差成比例减少”的形式。

这是比较模型。例如，为了比较从无条件模型到条件 1 模型在级别 1 中解释的方差成比例减少：

(intervals(unconditional.identity)$sigma[2]^2 - intervals(conditional1.identity)$sigma[2]^2) / intervals(unconditional.identity)$sigma[2]^2

希望您可以为您拥有的 2 级错误术语的数量“即插即用”相同的代码（在某些情况下不止一个）。确保以这种方式仅比较嵌套模型。

于 2017-12-08T23:00:52.577 回答