从SAS转换重复测量混合模型公式到R

Question

有一些关于混合模型的问题和帖子涉及更复杂的实验设计，因此我认为这个更简单的模型不仅可以帮助其他初学者，也可以帮助我自己。

所以，我的问题是我想从SAS PROC MIXED过程中制定一个R中的重复测量ancova：

proc mixed data=df1;
FitStatistics=akaike
class GROUP person day;
model Y = GROUP X1 / solution alpha=.1 cl;
repeated / type=cs subject=person group=GROUP;
lsmeans GROUP;
run;

这里是使用 R 创建的数据，输出的 SAS 结果如下所示：

.           Effect       panel    Estimate       Error      DF    t Value    Pr > |t|     Alpha       Lower       Upper
            Intercept              -9.8693      251.04       7      -0.04      0.9697       0.1     -485.49      465.75
            panel        1         -247.17      112.86       7      -2.19      0.0647       0.1     -460.99    -33.3510
            panel        2               0           .       .        .         .             .           .           .
            X1                     20.4125     10.0228       7       2.04      0.0811       0.1      1.4235     39.4016

以下是我使用'nlme'包在R中构建模型的方法，但是我没有得到类似的系数估计：

## create reproducible example fake panel data set:
set.seed(94); subject.id = abs(round(rnorm(10)*10000,0))

set.seed(99); sds = rnorm(10,15,5);means = 1:10*runif(10,7,13);trends = runif(10,0.5,2.5)

this = NULL; set.seed(98)
for(i in 1:10) { this = c(this,rnorm(6, mean = means[i], sd = sds[i])*trends[i]*1:6)}
set.seed(97)
that = sort(rep(rnorm(10,mean = 20, sd = 3),6))

df1 = data.frame(day = rep(1:6,10), GROUP = c(rep('TEST',30),rep('CONTROL',30)),
                 Y = this,
                 X1 = that,
                 person = sort(rep(subject.id,6)))

## use package nlme
require(nlme)

## run repeated measures mixed model using compound symmetry covariance structure:
summary(lme(Y ~ GROUP + X1, random = ~ +1 | person,
            correlation=corCompSymm(form=~day|person), na.action = na.exclude,
            data = df1,method='REML'))

现在，我认识到 R 的输出与 lm() 的输出类似：

                Value Std.Error DF    t-value p-value
(Intercept) -626.1622  527.9890 50 -1.1859379  0.2413
GROUPTEST   -101.3647  156.2940  7 -0.6485518  0.5373
X1            47.0919   22.6698  7  2.0772934  0.0764

我相信我已经接近规格了，但不确定我缺少什么东西使结果匹配（在合理范围内..）。任何帮助将不胜感激！

更新：使用下面答案中的代码，R输出变为：

> summary(model2)

滚动到底部查看参数估计值--看！与SAS完全相同。

Linear mixed-effects model fit by REML
 Data: df1 
      AIC      BIC   logLik
  776.942 793.2864 -380.471

Random effects:
 Formula: ~GROUP - 1 | person
 Structure: Diagonal
        GROUPCONTROL GROUPTEST Residual
StdDev:      184.692  14.56864 93.28885

Correlation Structure: Compound symmetry
 Formula: ~day | person 
 Parameter estimate(s):
         Rho 
-0.009929987 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | GROUP 
 Parameter estimates:
    TEST  CONTROL 
1.000000 3.068837

Fixed effects: Y ~ GROUP + X1 

                Value Std.Error DF    t-value p-value
(Intercept)   -9.8706 251.04678 50 -0.0393178  0.9688
GROUPTEST   -247.1712 112.85945  7 -2.1900795  0.0647
X1            20.4126  10.02292  7  2.0365914  0.0811

- baha-kev

---

你说的“没有得到相似的结果”是什么意思？是指有信息缺失，还是指得到了不同的估计值？如果是后者，你确定输入数据是一样的吗？ - David Robinson

这可能只是固定效应对比的差异吗？即 contrasts(df1$GROUP) <- contr.SAS(2)? - Ben Bolker

嗨@BenBolker！很高兴看到您注意到了这个帖子。我很好奇您是否同意我的下面的评估。我认为这比原帖作者期望的要棘手，但如果我错了那就太好了。 - Aaron left Stack Overflow

如果您添加统计信息并更多地询问哪种模型是适当的，那么这将是一个非常适合在stats.stackexchange上提问的问题。 - Aaron left Stack Overflow

嗨，亚伦，一旦我消化了这个问题，我会看看是否能以有趣且有教益的方式来表述问题，以引发该网站上的思考。 - baha-kev

Answer 1

请尝试以下操作：

model1 <- lme(
  Y ~ GROUP + X1,
  random = ~ GROUP | person,
  correlation = corCompSymm(form = ~ day | person),
  na.action = na.exclude, data = df1, method = "REML"
)
summary(model1)

我认为在这种情况下，使用 R 的 random = ~ groupvar | subjvar 选项与 SAS/MIXED 的 repeated / subject = subjvar group = groupvar 选项提供了类似的结果。

编辑:

SAS/MIXED

R (一个修订后的 model2)

model2 <- lme(
  Y ~ GROUP + X1,
  random = list(person = pdDiag(form = ~ GROUP - 1)),
  correlation = corCompSymm(form = ~ day | person),
  weights = varIdent(form = ~ 1 | GROUP),
  na.action = na.exclude, data = df1, method = "REML"
)
summary(model2)

因此，我认为这些协方差结构非常相似（σ_g1 = τ_g² + σ₁）。

编辑2：

协变量估计值（SAS/MIXED）：

Variance            person          GROUP TEST        8789.23
CS                  person          GROUP TEST         125.79
Variance            person          GROUP CONTROL       82775
CS                  person          GROUP CONTROL       33297

所以

TEST group diagonal element
  = 125.79 + 8789.23
  = 8915.02
CONTROL group diagonal element
  = 33297 + 82775
  = 116072

对角线元素 = σ_k1 + σ_k²。

协变量估计值（R lme）：

Random effects:
 Formula: ~GROUP - 1 | person
 Structure: Diagonal
        GROUP1TEST GROUP2CONTROL Residual
StdDev:   14.56864       184.692 93.28885

Correlation Structure: Compound symmetry
 Formula: ~day | person 
 Parameter estimate(s):
         Rho 
-0.009929987 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | GROUP 
 Parameter estimates:
   1TEST 2CONTROL 
1.000000 3.068837 

So

TEST group diagonal element
  = 14.56864^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + 93.28885^2
  = 8913.432
CONTROL group diagonal element
  = 184.692^2  + (3.068837^0.5 * 93.28885 * -0.009929987) + (3.068837 * 93.28885)^2
  = 116070.5

where diagonal element = τ_g² + σ₁ + σ_g².

- Triad sou.

---

我非常确定 random = ~ GROUP | person 不会改变原始代码，因为每个人只属于一个组。这种语法的作用是允许在个体内部组别水平之间的协方差有所不同。 - Aaron left Stack Overflow

我仍然认为随机效应不正确。random = list(person = pdDiag(form = ~ GROUP - 1)) 仍然允许组别水平之间的协方差在个体内不同，但是强制它们不相关。 - Aaron left Stack Overflow

另外，由于R实际上是通过相关性和方差来参数化模型的，因此很难看出您的矩阵如何与您编写的R代码相匹配。这并不一定是错误的，但如果您使用R的参数化来解释会更有帮助。 - Aaron left Stack Overflow

这几乎完美地匹配了proc mixed输出，所以我会标记为已回答。谢谢Triad！ - baha-kev

我只是尝试为 SAS/MIXED 中的 repeated / subject = subjvar group = groupvar 语句编写一个等效的 R 脚本。我同意相关线不是必要的。但是，在高度相关的数据中，没有相关线会导致失败（给出不同的结果）。可能需要进一步调查。 - Triad sou.

Answer 2

哦，这将是一个棘手的问题，如果要使用标准nlme函数实现这一点，需要认真研究Pinheiro/Bates的资料。但在你付出这些时间之前，你应该确保这是你所需的确切模型。也许有其他更适合您数据情况的方法。或者也许有些R可以更容易地完成同样好的工作，但不完全相同。首先，以下是我对您在SAS中进行的操作的理解：

repeated / type=cs subject=person group=GROUP;

这种type=cs subject=person会在同一个人的所有测量值之间引入相关性，而对于所有天数对，该相关性都是相同的。 group=GROUP 允许每个组的相关性不同。

相比之下，以下是我对您R代码正在执行的内容的看法：

random = ~ +1 | person,
correlation=corCompSymm(form=~day|person)

这段代码实际上以两种不同的方式添加了几乎相同的效果；random行为每个人添加了一个随机效果，而correlation行则在同一人的所有测量值之间引入了相关性。然而，这两件事情几乎是相同的；如果相关性是正的，那么包含其中任何一个都会得到完全相同的结果。我不确定如果同时包含两者会发生什么，但我知道只有一个是必要的。无论如何，此代码对所有个体具有相同的相关性，它不允许每个组拥有自己的相关性。

为了让每个组拥有自己的相关性，我认为您需要从两个不同的部分构建一个更复杂的相关结构；我从未做过这个，但我非常确定我记得Pinheiro/Bates曾经这样做过。

您可以考虑为每个人添加一个随机效应，然后使用weights=varIdent(form=~1|group)让不同组的方差不同（请检查我的语法是否正确）。这不完全相同，但讲述了类似的故事。SAS中的故事是某些个体的测量值更相关，而R中的故事是个体内的测量变异性不同；思考一下，具有更高变异性的测量将具有较低的相关性。因此，它们确实讲述了类似的故事，但是从相反的角度来看。

甚至可能（但我会感到惊讶）这两个模型最终成为相同事物的不同参数化。我的直觉是整体测量变异性在某种程度上会有所不同。但即使它们不是相同的东西，写出参数化也值得一试，以确保您理解它们并确保它们适当地描述了数据的情况。

- Aaron left Stack Overflow

---

最后一点，尽管如此，更改相关结构通常不会影响固定效应的估计值，因此如果有所不同，可能还有其他原因。 - Aaron left Stack Overflow

你的回答听起来对我来说很合理，但我真的认为我们需要从OP那里听/看到更多关于每个程序输出是什么样子的信息（我可以访问SAS，但不方便），以及主要差异在哪里... - Ben Bolker

感谢@BenBolker。我也没有尝试运行OP的代码； 我可以访问SAS，但在家中不方便。 - Aaron left Stack Overflow

一个想法--当我删除random=一词而保留correlation=一词时，我会遇到错误。如果像你说的那样，它们是多余的，那么我应该能够删除其中一个。 - baha-kev

@baha-kev：当您移除随机效应时，请确保切换到“gls”。 - Aaron left Stack Overflow