6

对于更复杂的实验设计，有几个关于混合模型的问题和帖子，所以我认为这个更简单的模型会帮助这个过程中的其他初学者以及我。

所以，我的问题是我想用 sas proc 混合程序在 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
4

2 回答 2

5

请尝试以下：

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选项 with提供了与选项R lme类似的结果。repeated / subject = subjvar group = groupvarSAS/MIXED

编辑：

SAS/混合

SAS/混合协方差矩阵

R（修改后的模型2）

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)

R协方差矩阵

所以，我认为这些协方差结构非常相似（σ g1 = τ g 2 + σ 1）。

编辑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 2。

协变量估计（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

其中对角元素 = τ g 2 + σ 1 + σ g 2。

于 2012-08-06T00:33:56.477 回答
4

哦，这将是一个棘手的问题，如果甚至可以使用标准的 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 中的故事是个体内部测量的可变性会有所不同。考虑到这一点，具有较高可变性的测量具有较低的相关性。所以他们确实讲述了类似的故事，但来自对立面。

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

于 2012-08-06T01:34:23.317 回答