将重复测量混合模型公式从 SAS 转换为 R [英] Converting Repeated Measures mixed model formula from SAS to R
问题描述
关于更复杂的实验设计的混合模型有几个问题和帖子，所以我认为这个更简单的模型可以帮助其他初学者以及我.
所以，我的问题是我想从 sas proc 混合程序中制定 R 中的重复测量 ancova:
proc 混合数据=df1;FitStatistics=akaike类 GROUP 人日；模型 Y = GROUP X1/解决方案 alpha=.1 cl；重复/类型=cs 主题=人组=组；lsmeans 组；跑步;这是使用在 R 中创建的数据的 SAS 输出(如下):
<预><代码>.效果面板估计误差 DF t 值 Pr >|t|Alpha 下 上截距 -9.8693 251.04 7 -0.04 0.9697 0.1 -485.49 465.75面板 1 -247.17 112.86 7 -2.19 0.0647 0.1 -460.99 -33.3510面板 2 0 .......X1 20.4125 10.0228 7 2.04 0.0811 0.1 1.4235 39.4016以下是我如何使用nlme"包在 R 中制定模型，但没有得到类似的系数估计:
## 创建可重现的示例假面板数据集:设置种子(94)；subject.id = abs(round(rnorm(10)*10000,0))设置种子(99)；sds = rnorm(10,15,5);means = 1:10*runif(10,7,13);trends = runif(10,0.5,2.5)这 = NULL;set.seed(98)for(i in 1:10) { this = c(this,rnorm(6, mean = mean[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 = 这个，X1 = 那，person = sort(rep(subject.id,6)))## 使用包 nlme要求(nlme)##使用复合对称协方差结构运行重复测量混合模型:总结(lme(Y ~ GROUP + X1，随机 = ~ +1 | 人，相关性=corCompSymm(form=~day|person), na.action = na.exclude,数据 = df1,method='REML'))现在，我现在意识到 R 的输出类似于 lm() 的输出:
Value Std.Error DF t-value p-value(拦截)-626.1622 527.9890 50 -1.1859379 0.2413组测试 -101.3647 156.2940 7 -0.6485518 0.5373X1 47.0919 22.6698 7 2.0772934 0.0764我相信我已经接近规范，但不确定我缺少哪一部分来使结果匹配(在合理范围内......).任何帮助将不胜感激！
更新: 使用下面答案中的代码，R 输出变为:
<代码>>总结(模型2)
滚动到底部的参数估计——看！与 SAS 相同.
REML 拟合的线性混合效应模型数据:df1AIC BIC logLik776.942 793.2864 -380.471随机效应:公式:~GROUP - 1 |人结构:对角线GROUPCONTROL GROUPTEST 残差标准差:184.692 14.56864 93.28885相关结构:复合对称公式:~日 |人参数估计:罗-0.009929987方差函数:结构:每层不同的标准差公式:~1 |团体参数估计:测控1.000000 3.068837固定效果:Y ~ GROUP + X1Value Std.Error DF t-value p-value(拦截) -9.8706 251.04678 50 -0.0393178 0.9688组测试 -247.1712 112.85945 7 -2.1900795 0.0647X1 20.4126 10.02292 7 2.0365914 0.0811请尝试以下:
model1 <- lme(Y~组+X1，随机 = ~ 组 |人，相关性 = corCompSymm(形式 = ~ 天 | 人)，na.action = na.exclude, 数据 = df1, 方法 = "REML")总结(模型1)我认为 random = ~ groupvar |subjvar 选项和 R lme 提供了与 repeated/subject = subjvar group = groupvar 选项和 SAS/MIXED 类似的结果 在这种情况下.
SAS/混合
R(修改后的模型 2)
model2 <- lme(Y~组+X1，random = list(person = pdDiag(form = ~ GROUP - 1)),相关性 = corCompSymm(形式 = ~ 天 | 人)，weights = varIdent(form = ~ 1 | GROUP),na.action = na.exclude, 数据 = df1, 方法 = "REML")总结(模型2)所以，我认为这些协方差结构非常相似(σg1 = τg2 + σ1).
编辑 2:
协变量估计值 (SAS/MIXED):
方差人组测试 8789.23CS人组测试125.79差异人 GROUP CONTROL 82775CS 人 群控 33297所以
TEST组对角线元素= 125.79 + 8789.23= 8915.02CONTROL 组对角线元素= 33297 + 82775= 116072其中对角线元素 = σk1 + σk2.
协变量估计 (R lme):
随机效果:公式:~GROUP - 1 |人结构:对角线GROUP1TEST GROUP2CONTROL 残差标准差:14.56864 184.692 93.28885相关结构:复合对称公式:~日 |人参数估计:罗-0.009929987方差函数:结构:每层不同的标准差公式:~1 |团体参数估计:1测试 2控制1.000000 3.068837所以
TEST组对角线元素= 14.56864^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + 93.28885^2= 8913.432CONTROL 组对角线元素= 184.692^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + (3.068837 * 93.28885)^2= 116070.5其中对角线元素 = τg2 + σ1 + σg2.
There are several questions and posts about mixed models for more complex experimental designs, so I thought this more simple model would help other beginners in this process as well as I.
So, my question is I would like to formulate a repeated measures ancova in R from sas proc mixed procedure:
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;
Here is the SAS output using the data created in R (below):
. 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
Below is how I formulated the model in R using 'nlme' package, but am not getting similar coefficient estimates:
## 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'))
Now, the output from R, which I now realize is similar to the output from 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
I believe I'm close as to the specification, but not sure what piece I'm missing to make the results match (within reason..). Any help would be appreciated!
UPDATE: Using the code in the answer below, the R output becomes:
> summary(model2)
Scroll to bottom for the parameter estimates -- look! identical to 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
Please try below:
model1 <- lme(
Y ~ GROUP + X1,
random = ~ GROUP | person,
correlation = corCompSymm(form = ~ day | person),
na.action = na.exclude, data = df1, method = "REML"
)
summary(model1)
I think random = ~ groupvar | subjvar option with R lme provides similar result of repeated / subject = subjvar group = groupvar option with SAS/MIXED in this case.
Edit:
SAS/MIXED
R (a revised 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)
So, I think these covariance structures are very similar (σg1 = τg2 + σ1).
Edit 2:
Covariate estimates (SAS/MIXED):
Variance person GROUP TEST 8789.23
CS person GROUP TEST 125.79
Variance person GROUP CONTROL 82775
CS person GROUP CONTROL 33297
So
TEST group diagonal element
= 125.79 + 8789.23
= 8915.02
CONTROL group diagonal element
= 33297 + 82775
= 116072
where diagonal element = σk1 + σk2.
Covariate estimates (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 = τg2 + σ1 + σg2.
这篇关于将重复测量混合模型公式从 SAS 转换为 R的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
