获取具有随机截距和具有异质方差的连续一阶自回归相关结构的 varcov 矩阵 LMM？
Get varcov matrix LMM with random intercept and continuous first order autoregressive correlation structure with heteregenous variances?
我有 66 名内源性或外源性抑郁症患者 (endo) 的重复测量数据和每周测量的抑郁评分，持续 0-5 周 (hdrs 和 week，所以每个患者六次测量，包括基线）。数据为长格式：
library(dplyr)
library(magrittr)
library(nlme)
mydata <- structure(list(id = c(101, 101, 101, 101, 101, 101, 103, 103,
103, 103, 103, 103, 104, 104, 104, 104, 104, 104, 105, 105, 105,
105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107,
108, 108, 108, 108, 108, 108, 113, 113, 113, 113, 113, 114, 114,
114, 114, 114, 115, 115, 115, 115, 115, 117, 117, 117, 117, 117,
117, 118, 118, 118, 118, 118, 120, 120, 120, 120, 120, 120, 121,
121, 121, 121, 121, 121, 123, 123, 123, 123, 123, 123, 302, 302,
302, 302, 302, 302, 303, 303, 303, 303, 303, 303, 304, 304, 304,
304, 304, 305, 305, 305, 305, 305, 305, 308, 308, 308, 308, 308,
308, 309, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311,
311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313,
313, 313, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 316,
318, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 319, 322,
322, 322, 322, 322, 327, 327, 327, 327, 327, 327, 328, 328, 328,
328, 328, 328, 331, 331, 331, 331, 331, 331, 333, 333, 333, 333,
333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 335,
337, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 338, 339,
339, 339, 339, 339, 344, 344, 344, 344, 344, 345, 345, 345, 345,
345, 345, 346, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347,
348, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 349, 350,
350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 351, 352, 352,
352, 352, 352, 352, 353, 353, 353, 353, 353, 353, 354, 354, 354,
354, 354, 355, 355, 355, 355, 355, 355, 357, 357, 357, 357, 357,
357, 360, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 361,
501, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 502, 504,
504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 505, 507, 507,
507, 507, 507, 507, 603, 603, 603, 603, 603, 603, 604, 604, 604,
604, 604, 606, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607,
607, 608, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610,
610, 610, 610), week = structure(c(0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0,
1, 2, 3, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 1, 2, 3, 4, 5, 1,
2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5,
0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 5, 0, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4,
5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 2, 3, 4,
5, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3,
4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1,
2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4,
5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1,
2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5,
1, 2, 3, 4, 5, 0, 2, 3, 5), format.spss = "F1.0", display_width = 6L),
week_fact = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 6L,
1L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L,
5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 3L, 4L, 5L,
6L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L,
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L,
5L, 6L, 1L, 3L, 4L, 6L), .Label = c("Week 0", "Week 1", "Week 2",
"Week 3", "Week 4", "Week 5"), class = "factor"), endo = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Exogenous",
"Endogenous"), class = "factor"), hdrs = c(26, 22, 18, 7,
4, 3, 33, 24, 15, 24, 15, 13, 29, 22, 18, 13, 19, 0, 22,
12, 16, 16, 13, 9, 21, 25, 23, 18, 20, 21, 21, 16, 19, 6,
21, 22, 11, 9, 9, 7, 21, 23, 19, 23, 23, 17, 11, 13, 7, 7,
16, 16, 16, 16, 11, 19, 16, 13, 12, 7, 6, 26, 18, 18, 14,
11, 20, 19, 17, 18, 16, 17, 20, 22, 19, 19, 12, 14, 15, 15,
15, 13, 5, 5, 18, 22, 16, 8, 9, 12, 21, 21, 13, 14, 10, 5,
21, 27, 29, 12, 24, 19, 17, 15, 11, 5, 1, 22, 21, 18, 17,
12, 11, 22, 22, 16, 19, 20, 11, 24, 19, 11, 7, 6, 20, 16,
21, 17, 15, 17, 18, 17, 17, 6, 21, 19, 10, 11, 11, 8, 27,
21, 17, 13, 5, 32, 26, 23, 26, 23, 24, 17, 18, 19, 21, 17,
11, 24, 18, 10, 14, 13, 12, 28, 21, 25, 32, 34, 17, 18, 15,
8, 19, 17, 22, 24, 28, 26, 28, 29, 19, 21, 18, 16, 14, 10,
23, 20, 21, 20, 24, 14, 31, 25, 7, 8, 11, 21, 21, 18, 15,
12, 10, 27, 22, 23, 21, 12, 13, 22, 20, 22, 23, 19, 18, 27,
14, 12, 11, 12, 21, 12, 13, 13, 18, 29, 27, 27, 22, 22, 23,
25, 24, 19, 23, 14, 21, 18, 15, 14, 10, 8, 24, 21, 12, 13,
12, 5, 17, 19, 15, 12, 9, 13, 22, 25, 12, 16, 10, 16, 30,
27, 23, 20, 12, 11, 21, 19, 18, 15, 18, 19, 27, 21, 24, 22,
16, 11, 28, 27, 27, 26, 23, 22, 26, 20, 13, 10, 7, 27, 22,
24, 25, 19, 19, 21, 28, 27, 29, 28, 33, 30, 22, 11, 8, 7,
19, 29, 30, 26, 22, 19, 24, 21, 22, 13, 11, 2, 1, 19, 17,
15, 16, 12, 12, 21, 11, 18, 0, 0, 4, 27, 26, 26, 25, 24,
19, 28, 22, 18, 20, 11, 13, 27, 27, 13, 5, 7, 19, 33, 12,
12, 3, 1, 30, 39, 30, 27, 20, 4, 24, 19, 14, 12, 3, 4, 25,
22, 14, 15, 2, 34, 33, 23, 11)), row.names = c(NA, -375L), class = "data.frame")
head(mydata)
id week week_fact endo hdrs
1 101 0 Week 0 Exogenous 26
2 101 1 Week 1 Exogenous 22
3 101 2 Week 2 Exogenous 18
4 101 3 Week 3 Exogenous 7
5 101 4 Week 4 Exogenous 4
6 101 5 Week 5 Exogenous 3
为了检查抑郁症的变化模式 (hdrs) 在各组之间是否不同，我 运行 有两个模型。第一个包括时间的随机截距和斜率。第二个模型包括时间的随机截距（无随机斜率）并假设 连续一阶自回归相关结构，残差具有异质方差 :
model_ris <-
lme(fixed=hdrs ~ endo*week,
random=~1 + week|id,
method="ML",
data=mydata)
model_ri_car1_hetvar <-
gls(hdrs ~ endo*week,
correlation=corCAR1(form=~1|id),
weight=varIdent(form=~1|week),
method="ML",
data=mydata)
我有两个问题：
我是否正确指定了第二个模型？我从两个不同的模型中借用了代码（correlation=corCAR1(form=~1|id) 部分来自使用 gls 指定 cCAR(1) 模型的人，weight=varIdent(form=~1|week) 来自 lme 最初在具有异质方差模型的 AR(1) 中），但我非常不确定这是否是我指定我正在尝试做的事情的方式。
有没有办法得到第二个模型的方差协方差矩阵？我试过 getVarCov 但这不起作用。调用 corMatrix(model_ri_car1_hetvar) 给出了一个相关矩阵，但我想要一个方差-协方差矩阵。
如有任何帮助，我们将不胜感激！
谢谢！
我不确定您为什么从 lme 模型切换到 gls 模型。你的描述说“[t]他的第二个模型包括时间的随机拦截”，但是（a）我什至不确定那是什么意思（你的意思是说“对患者的随机拦截（id） " ?) 和 (b) gls 模型 不实现随机效应 （它们用于当你想要建模相关性和异方差并且 不't 想要在模型中包含随机效应）。我会说如果你想要一个随机截距（跨患者），一个连续时间自相关结构（在患者内），以及每周不同的方差，你应该使用
model_ri_car1_hetvar <-
update(model_ris,
random = ~1 |id,
correlation=corCAR1(form=~1|id),
weight=varIdent(form=~1|week))
（仅使用已更改的组件更新模型）。
为了得到协方差矩阵，我会说你应该：（1）select每周测量的某个人的相关矩阵；
cs <- model_ri_car1_hetvar$modelStruct$corStruct
cormat <- corMatrix(cs)[["101"]] ## this individual happens be measured every week
(2) 得到每周的标准差：
hs <- model_ri_car1_hetvar$modelStruct$varStruct
## ?varIdent: "the coefficients of the variance function represent
## the ratios between the variances and a reference variance
## (corresponding to a reference group level).
sdvec <- sigma(model_ri_car1_hetvar)*c(1, sqrt(coef(hs, unconstrained = FALSE)))
(3) 计算协方差矩阵：
diag(sdvec) %*% cormat %*% diag(sdvec)
在一些模拟数据上仔细检查everything/try这不会有什么坏处！
来自你的简历问题：
How much parameters for the correlation structure of residuals does the second model actually estimate? I know that AR(1) and cAR(1) models just estimate phi(or rho), but now that I've combined these correlation structures with heterogeneous variances, I'm not sure if my model is estimating a rho/phi for every timepoint for every patient (which wouldn't be a very good option). Could I see/check this anywhere?
一个 phi 值。当您 运行 summary() 时，您可以看到参数估计值 (Phi: 0.5808927)。或者，coef(cs, unconstrained = FALSE) 为您提供该值。
当我向模型添加随机斜率项时：
model_ris_car1_hetvar <-
update(model_ri_car1_hetvar,
random = ~1 + week |id)
它确实有效（R-devel，Linux）：“奇异收敛”的东西往往是非常特定于平台和版本的。我会尝试 control = lmeControl(opt = "optim") 看看是否有帮助...
比较模型，随机斜率似乎值得保留，但相关性和异方差性可能不值得（取决于您的建模目标）。
bbmle::AICtab(model_ris, model_ri_car1_hetvar, model_ris_car1_hetvar)
dAIC df
model_ris 0.0 8
model_ris_car1_hetvar 0.5 14
model_ri_car1_hetvar 11.4 12
我有 66 名内源性或外源性抑郁症患者 (endo) 的重复测量数据和每周测量的抑郁评分，持续 0-5 周 (hdrs 和 week，所以每个患者六次测量，包括基线）。数据为长格式：
library(dplyr)
library(magrittr)
library(nlme)
mydata <- structure(list(id = c(101, 101, 101, 101, 101, 101, 103, 103,
103, 103, 103, 103, 104, 104, 104, 104, 104, 104, 105, 105, 105,
105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107,
108, 108, 108, 108, 108, 108, 113, 113, 113, 113, 113, 114, 114,
114, 114, 114, 115, 115, 115, 115, 115, 117, 117, 117, 117, 117,
117, 118, 118, 118, 118, 118, 120, 120, 120, 120, 120, 120, 121,
121, 121, 121, 121, 121, 123, 123, 123, 123, 123, 123, 302, 302,
302, 302, 302, 302, 303, 303, 303, 303, 303, 303, 304, 304, 304,
304, 304, 305, 305, 305, 305, 305, 305, 308, 308, 308, 308, 308,
308, 309, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311,
311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313,
313, 313, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 316,
318, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 319, 322,
322, 322, 322, 322, 327, 327, 327, 327, 327, 327, 328, 328, 328,
328, 328, 328, 331, 331, 331, 331, 331, 331, 333, 333, 333, 333,
333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 335,
337, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 338, 339,
339, 339, 339, 339, 344, 344, 344, 344, 344, 345, 345, 345, 345,
345, 345, 346, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347,
348, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 349, 350,
350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 351, 352, 352,
352, 352, 352, 352, 353, 353, 353, 353, 353, 353, 354, 354, 354,
354, 354, 355, 355, 355, 355, 355, 355, 357, 357, 357, 357, 357,
357, 360, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 361,
501, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 502, 504,
504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 505, 507, 507,
507, 507, 507, 507, 603, 603, 603, 603, 603, 603, 604, 604, 604,
604, 604, 606, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607,
607, 608, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610,
610, 610, 610), week = structure(c(0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0,
1, 2, 3, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 1, 2, 3, 4, 5, 1,
2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5,
0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 5, 0, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4,
5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 2, 3, 4,
5, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2,
3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0,
1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3,
4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1,
2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4,
5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1,
2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5,
1, 2, 3, 4, 5, 0, 2, 3, 5), format.spss = "F1.0", display_width = 6L),
week_fact = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 6L,
1L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L,
5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 3L, 4L, 5L,
6L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L,
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L,
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L,
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L,
6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L,
5L, 6L, 1L, 3L, 4L, 6L), .Label = c("Week 0", "Week 1", "Week 2",
"Week 3", "Week 4", "Week 5"), class = "factor"), endo = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Exogenous",
"Endogenous"), class = "factor"), hdrs = c(26, 22, 18, 7,
4, 3, 33, 24, 15, 24, 15, 13, 29, 22, 18, 13, 19, 0, 22,
12, 16, 16, 13, 9, 21, 25, 23, 18, 20, 21, 21, 16, 19, 6,
21, 22, 11, 9, 9, 7, 21, 23, 19, 23, 23, 17, 11, 13, 7, 7,
16, 16, 16, 16, 11, 19, 16, 13, 12, 7, 6, 26, 18, 18, 14,
11, 20, 19, 17, 18, 16, 17, 20, 22, 19, 19, 12, 14, 15, 15,
15, 13, 5, 5, 18, 22, 16, 8, 9, 12, 21, 21, 13, 14, 10, 5,
21, 27, 29, 12, 24, 19, 17, 15, 11, 5, 1, 22, 21, 18, 17,
12, 11, 22, 22, 16, 19, 20, 11, 24, 19, 11, 7, 6, 20, 16,
21, 17, 15, 17, 18, 17, 17, 6, 21, 19, 10, 11, 11, 8, 27,
21, 17, 13, 5, 32, 26, 23, 26, 23, 24, 17, 18, 19, 21, 17,
11, 24, 18, 10, 14, 13, 12, 28, 21, 25, 32, 34, 17, 18, 15,
8, 19, 17, 22, 24, 28, 26, 28, 29, 19, 21, 18, 16, 14, 10,
23, 20, 21, 20, 24, 14, 31, 25, 7, 8, 11, 21, 21, 18, 15,
12, 10, 27, 22, 23, 21, 12, 13, 22, 20, 22, 23, 19, 18, 27,
14, 12, 11, 12, 21, 12, 13, 13, 18, 29, 27, 27, 22, 22, 23,
25, 24, 19, 23, 14, 21, 18, 15, 14, 10, 8, 24, 21, 12, 13,
12, 5, 17, 19, 15, 12, 9, 13, 22, 25, 12, 16, 10, 16, 30,
27, 23, 20, 12, 11, 21, 19, 18, 15, 18, 19, 27, 21, 24, 22,
16, 11, 28, 27, 27, 26, 23, 22, 26, 20, 13, 10, 7, 27, 22,
24, 25, 19, 19, 21, 28, 27, 29, 28, 33, 30, 22, 11, 8, 7,
19, 29, 30, 26, 22, 19, 24, 21, 22, 13, 11, 2, 1, 19, 17,
15, 16, 12, 12, 21, 11, 18, 0, 0, 4, 27, 26, 26, 25, 24,
19, 28, 22, 18, 20, 11, 13, 27, 27, 13, 5, 7, 19, 33, 12,
12, 3, 1, 30, 39, 30, 27, 20, 4, 24, 19, 14, 12, 3, 4, 25,
22, 14, 15, 2, 34, 33, 23, 11)), row.names = c(NA, -375L), class = "data.frame")
head(mydata)
id week week_fact endo hdrs
1 101 0 Week 0 Exogenous 26
2 101 1 Week 1 Exogenous 22
3 101 2 Week 2 Exogenous 18
4 101 3 Week 3 Exogenous 7
5 101 4 Week 4 Exogenous 4
6 101 5 Week 5 Exogenous 3
为了检查抑郁症的变化模式 (hdrs) 在各组之间是否不同，我 运行 有两个模型。第一个包括时间的随机截距和斜率。第二个模型包括时间的随机截距（无随机斜率）并假设 连续一阶自回归相关结构，残差具有异质方差 :
model_ris <-
lme(fixed=hdrs ~ endo*week,
random=~1 + week|id,
method="ML",
data=mydata)
model_ri_car1_hetvar <-
gls(hdrs ~ endo*week,
correlation=corCAR1(form=~1|id),
weight=varIdent(form=~1|week),
method="ML",
data=mydata)
我有两个问题：
我是否正确指定了第二个模型？我从两个不同的模型中借用了代码（
correlation=corCAR1(form=~1|id)部分来自使用gls指定 cCAR(1) 模型的人，weight=varIdent(form=~1|week)来自lme最初在具有异质方差模型的 AR(1) 中），但我非常不确定这是否是我指定我正在尝试做的事情的方式。有没有办法得到第二个模型的方差协方差矩阵？我试过
getVarCov但这不起作用。调用corMatrix(model_ri_car1_hetvar)给出了一个相关矩阵，但我想要一个方差-协方差矩阵。
如有任何帮助，我们将不胜感激！ 谢谢！
我不确定您为什么从 lme 模型切换到 gls 模型。你的描述说“[t]他的第二个模型包括时间的随机拦截”，但是（a）我什至不确定那是什么意思（你的意思是说“对患者的随机拦截（id） " ?) 和 (b) gls 模型 不实现随机效应 （它们用于当你想要建模相关性和异方差并且 不't 想要在模型中包含随机效应）。我会说如果你想要一个随机截距（跨患者），一个连续时间自相关结构（在患者内），以及每周不同的方差，你应该使用
model_ri_car1_hetvar <-
update(model_ris,
random = ~1 |id,
correlation=corCAR1(form=~1|id),
weight=varIdent(form=~1|week))
（仅使用已更改的组件更新模型）。
为了得到协方差矩阵，我会说你应该：（1）select每周测量的某个人的相关矩阵；
cs <- model_ri_car1_hetvar$modelStruct$corStruct
cormat <- corMatrix(cs)[["101"]] ## this individual happens be measured every week
(2) 得到每周的标准差：
hs <- model_ri_car1_hetvar$modelStruct$varStruct
## ?varIdent: "the coefficients of the variance function represent
## the ratios between the variances and a reference variance
## (corresponding to a reference group level).
sdvec <- sigma(model_ri_car1_hetvar)*c(1, sqrt(coef(hs, unconstrained = FALSE)))
(3) 计算协方差矩阵：
diag(sdvec) %*% cormat %*% diag(sdvec)
在一些模拟数据上仔细检查everything/try这不会有什么坏处！
来自你的简历问题：
How much parameters for the correlation structure of residuals does the second model actually estimate? I know that AR(1) and cAR(1) models just estimate phi(or rho), but now that I've combined these correlation structures with heterogeneous variances, I'm not sure if my model is estimating a rho/phi for every timepoint for every patient (which wouldn't be a very good option). Could I see/check this anywhere?
一个 phi 值。当您 运行 summary() 时，您可以看到参数估计值 (Phi: 0.5808927)。或者，coef(cs, unconstrained = FALSE) 为您提供该值。
当我向模型添加随机斜率项时：
model_ris_car1_hetvar <-
update(model_ri_car1_hetvar,
random = ~1 + week |id)
它确实有效（R-devel，Linux）：“奇异收敛”的东西往往是非常特定于平台和版本的。我会尝试 control = lmeControl(opt = "optim") 看看是否有帮助...
比较模型，随机斜率似乎值得保留，但相关性和异方差性可能不值得（取决于您的建模目标）。
bbmle::AICtab(model_ris, model_ri_car1_hetvar, model_ris_car1_hetvar)
dAIC df
model_ris 0.0 8
model_ris_car1_hetvar 0.5 14
model_ri_car1_hetvar 11.4 12