转换混合模型公式从SAS至R

Question

重复测量有几个问题和有关混合模式为更复杂的实验设计的职位，所以我想这更简单的模型将帮助其他新手则在这个过程中，以及一

所以，我的问题是我想从SAS制定中的R重复测量的协方差分析PROC MIXED程序：

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 

下面是我如何制定的R模型中使用“NLME”包，但我没有得到相似系数估计：

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

我相信我接近的规范，但不知道我错过了什么片，使结果匹配（在合理范围内。）。任何帮助，将不胜感激！

UPDATE：在下面的答案使用代码，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 

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) 

我觉得random = ~ groupvar | subjvar选项与Rlme提供了在这种情况下repeated/subject = subjvar group = groupvar选项与SAS/MIXED类似的结果。

编辑：

SAS /混合

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) 

所以，我觉得这些COVAR iance结构非常相似（σ _G1 τ =_克 + σ ）。

编辑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 + σ _ķ 。

协方差估计（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 

所以

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 

其中对角元素= 克 + σ + 克 。

Answer 2

哦，这将是一个棘手的问题，如果甚至有可能使用标准的nlme函数，将会对Pinheiro/Bates进行一些认真的研究。

在你花时间做这件事之前，你应该确保这是你需要的确切模型。也许还有别的东西可以更好地适合你的数据的故事。或者，也许R可以做得更容易，这是一样的好，但不完全一样。

首先，这里就你在SAS做这一行我的看法：

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

这type=cs subject=person是对同一人所有的测量值之间引发的相关性，而且相关性都是一样的双日。 允许每个组的相关性不同。

相反，这是我拿上你的R代码里面做：

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

此代码实际上是将几乎以两种不同的方式同样的效果; random行会为每个人添加一个随机效果，并且correlation行会导致同一个人的所有测量值之间的相关性。但是，这两件事情几乎完全相同。如果相关性为正数，则通过包含其中任何一个来得到完全相同的结果。我不确定当你同时包括这两个时会发生什么，但我知道只有一个是必要的。无论如何，这些代码对于所有个体都具有相同的相关性，它不允许每个群体拥有自己的相关性。

为了让每个小组都有自己的关联，我认为你必须在两个不同的部分之间建立一个更复杂的关联结构;我从来没有这样做过，但我很确定我记得Pinheiro/Bates在做这件事。

您可能会考虑为个人添加一个随机效果，然后让weights=varIdent(form=~1|group)（来自内存，请检查我的语法）的不同组的差异有所不同。这不完全相同，但讲述了类似的故事。 SAS中的故事是，对某些个体的测量比对其他个体的测量更加相关。思考这意味着什么，对于具有较高相关性的个体的测量结果将比对具有较低相关性的个体的测量结果更接近。相反，R中的故事是个体内测量的变化性不同;考虑到这一点，具有较高变异性的测量具有较低的相关性。所以他们确实讲述类似的故事，但从相反的方面来看。

这甚至有可能（但我会感到惊讶的），这两个模型最终成为同一事物的不同参数化。我的直觉是总体测量的可变性在某些方面会有所不同。但即使它们不是同一件事，为了确保您理解它们并确保它们适当地描述数据的故事，值得写出参数化值。