I have repeated measurements data on 66 patients with either endogenous or exogenous depression (endo) and depression scores measured weekly for 0-5 weeks (hdrs and week, so six measurements per patients including baseline). The data is in a long format:
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
To examine if the pattern of change of depression (hdrs) is different between groups I've run two models. The first includes a random intercept and slope for time. The second model includes a random intercept for time (no random slope) and assumes a continuous first order autoregressive correlation structure with heterogeneous variances for the residuals:
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)
I've got two questions:
Have I specified the second model correctly? I've borrowed code from two different models (the correlation=corCAR1(form=~1|id) part is from someone useing gls to specify a cCAR(1) model, and the weight=varIdent(form=~1|week) is from code from lme originally in a AR(1) with heterogenous variances model) but I'm very much not sure if this is how I specify what I'm trying to do.
Is there a way to get a variance covariance matrix for the second model? I've tried getVarCov but that does not work. Calling corMatrix(model_ri_car1_hetvar) gives a correlation matrix, but I want a variance-covariance matrix.
Any help would be greatly appreciated! Thanks!
I'm not sure why you've switched from an lme model to a gls model. Your description says "[t]he second model includes a random intercept for time", but (a) I'm not even sure what that means (did you mean to say "a random intercept across patients (id)" ?) and (b) gls models don't implement random effects (they're intended for when you want to model correlation and heteroscedasticity and don't want to include random effects in the model). I would say if you want a random intercept (across patients), a continuous-time autocorrelation structure (within patients), and different variances by week, you should use
model_ri_car1_hetvar <-
update(model_ris,
random = ~1 |id,
correlation=corCAR1(form=~1|id),
weight=varIdent(form=~1|week))
(updating model only with the components that have changed).
In order to get a covariance matrix, I would say you should: (1) select the correlation matrix for some individual who was measured in every week;
cs <- model_ri_car1_hetvar$modelStruct$corStruct
cormat <- corMatrix(cs)[["101"]] ## this individual happens be measured every week
(2) get the standard deviations for each week:
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) compute the covariance matrix:
diag(sdvec) %*% cormat %*% diag(sdvec)
It wouldn't hurt to double-check everything/try this on some simulated data!
From your CV question:
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?
One phi value. When you run summary() you can see the parameter estimates (Phi: 0.5808927). Alternately, coef(cs, unconstrained = FALSE) gives you that value.
When I add a random slope term to the model:
model_ris_car1_hetvar <-
update(model_ri_car1_hetvar,
random = ~1 + week |id)
it actually works (R-devel, Linux): the "singular convergence" stuff tends to be very platform- and version-specific. I would try control = lmeControl(opt = "optim") and see if it helps ...
Comparing the models, it seems that the random slopes are worth keeping but the correlations and heteroscedasticity might not be (depending on your modeling goal).
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