R与SPSS混合模型重复测量代码[来自交叉验证]

注意：这个问题最初发布在交叉验证中，建议在stackoverflow中询问它。

我正在尝试建模三路重复测量实验，fixedfactora * fixedFactorb *时间[天]。没有丢失的观察结果，但是我的小组（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中尝试的，它产生的AIC为2156：

    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)

这些是我的问题：

上面的SPSS代码产生了一个具有AIC = 2471的模型，而R代码产生了一个具有AIC = 2156的模型。使模型不同的代码是什么？
从我上面描述的内容中，是否适合我要测试的模型？如果没有，什么是更好的方法，我该如何在两个程序中获得相同的结果？

编辑

要注意的另一件事是，我没有虚拟编码我的因素。我不知道这是软件的问题，还是内置的编码在SPSS vs R中是否有所不同。我也不知道这是否是我三向互动术语的问题。p>另外，当我说"因素"时，我的意思是一个不变的群体或特征（例如"性"）。

从一个无条件模型开始，一个模型在1级时具有身份方差 - 稳定性结构，一个具有AR（1）var-covar结构在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')

找到该无条件模型的类内相关系数，该模型是2级误差除以级别1和2级误差的总和。在电子表格程序中，这可能更容易，但是在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)

这取决于您的领域，但是在教育研究中，ICC低于0.2，绝对低于0.1，被认为尚未准备好用于层次线性模型。也就是说，由于确认了独立性的假设，因此多重回归会更好。如果您的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')

如果您的因素是时间不变的（您在交叉验证中说），则您的模型会更大，因为时间和组都嵌套在这些固定效果中：

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()的系数上获得置信区间，或使用summary()的P值获得置信区。请记住，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级错误项的数量"插入并播放"相同的代码（在某些情况下是不止一个）。确保仅以这种方式比较嵌套模型。

相关内容

最新更新

热门标签：