I´m trying to fit longitudinal time to event data with a joint model using the R-package "JM". This is my first attempt with joint models and follows a textbook approach:
aids.id <- aids[!duplicated(aids$patient),]
lmeFit.aids <- lme(CD4~obstime + obstime:drug, random=~obstime|patient, data=aids)
coxFit.aids <- coxph(Sdurv(Time,death)~drug,data=aids.id, x=TRUE)
jointFit.aids <- jointModel(lmeFit.aids, coxFit.aids, timeVar="obstime",method="piecewise-PH-aGH")
summary(jointFit.aids)
The code works as expected. But when I use my own data, it doesn't work anymore.
str(DATA) 'data.frame': 6436 obs. of 13 variables: $ patnr : Factor w/ 1669 levels "0010000158","0010000278",..: 4 4 4 4 7 7 7 7 7 7 ...
$ sex : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ date : POSIXct, format: "2008-08-08" "2010-01-25" "2012-02-24" "2012-04-21" ...
$ time : num 1355 1355 1355 1355 834 ...
$ Crea : num 7.4 9.6 12.3 10.3 0.8 ...
$ CysC : num 6.2 5.84 6.17 5.32 0.9 0.94 0.92 0.91 0.91 0.91 ...
$ deathdate : POSIXct, format: "2012-04-24" "2012-04-24" "2012-04-24" "2012-04-24" ...
$ start_date: POSIXct, format: "2008-08-08" "2008-08-08" "2008-08-08" "2008-08-08" ...
$ stop_date : POSIXct, format: "2010-01-25" "2012-02-24" "2012-04-21" "2012-04-24" ...
$ start : num 0 535 1295 1352 0 ...
$ stop : num 535 1295 1352 1355 3 ...
$ obstime : num 0 535 1295 1352 0 ...
$ event : num 0 0 0 1 0 0 0 0 0 0 ...
These are the first 20 lines of the data frame:
patnr sex date time Crea CysC deathdate start_date stop_date start stop obstime event
637 0010000343 1 2008-08-08 1355.00000 7.40 6.20 2012-04-24 2008-08-08 2010-01-25 0.00000 535.04167 0.00000 0
816 0010000343 1 2010-01-25 1355.00000 9.60 5.84 2012-04-24 2008-08-08 2012-02-24 535.04167 1295.04167 535.04167 0
1171 0010000343 1 2012-02-24 1355.00000 12.31 6.17 2012-04-24 2008-08-08 2012-04-21 1295.04167 1352.00000 1295.04167 0
1201 0010000343 1 2012-04-21 1355.00000 10.35 5.32 2012-04-24 2008-08-08 2012-04-24 1352.00000 1355.00000 1352.00000 1
1363 0010000873 1 2011-12-05 834.00000 0.80 0.90 <NA> 2011-12-05 2011-12-08 0.00000 3.00000 0.00000 0
1370 0010000873 1 2011-12-08 834.00000 0.52 0.94 <NA> 2011-12-05 2011-12-09 3.00000 4.00000 3.00000 0
1372 0010000873 1 2011-12-09 834.00000 0.45 0.92 <NA> 2011-12-05 2011-12-18 4.00000 13.00000 4.00000 0
1386 0010000873 1 2011-12-18 834.00000 0.34 0.91 <NA> 2011-12-05 2011-12-19 13.00000 14.00000 13.00000 0
1387 0010000873 1 2011-12-19 834.00000 0.31 0.91 <NA> 2011-12-05 2011-12-20 14.00000 15.00000 14.00000 0
1391 0010000873 1 2011-12-20 834.00000 0.62 0.91 <NA> 2011-12-05 2011-12-27 15.00000 22.00000 15.00000 0
1411 0010000873 1 2011-12-27 834.00000 0.61 1.44 <NA> 2011-12-05 2011-12-31 22.00000 26.00000 22.00000 0
1418 0010000873 1 2011-12-31 834.00000 0.43 1.18 <NA> 2011-12-05 2012-01-01 26.00000 27.00000 26.00000 0
1419 0010000873 1 2012-01-01 834.00000 0.46 1.22 <NA> 2011-12-05 2013-07-09 27.00000 581.95833 27.00000 0
1466 0010000873 1 2013-07-09 834.00000 0.85 0.91 <NA> 2011-12-05 2014-03-18 581.95833 834.00000 581.95833 0
1478 0010000873 1 2014-03-18 834.00000 1.20 1.00 <NA> 2011-12-05 2015-09-18 834.00000 1382.95833 834.00000 0
2020 0010002412 1 2015-03-26 23.95833 1.16 0.85 <NA> 2015-03-26 2015-04-10 0.00000 14.95833 0.00000 0
2035 0010002412 1 2015-04-10 23.95833 0.67 0.74 <NA> 2015-03-26 2015-04-14 14.95833 18.95833 14.95833 0
2043 0010002412 1 2015-04-14 23.95833 0.56 0.75 <NA> 2015-03-26 2015-04-16 18.95833 20.95833 18.95833 0
2046 0010002412 1 2015-04-16 23.95833 0.45 0.75 <NA> 2015-03-26 2015-04-17 20.95833 21.95833 20.95833 0
2049 0010002412 1 2015-04-17 23.95833 0.52 0.86 <NA> 2015-03-26 2015-04-18 21.95833 22.95833 21.95833 0
This is the code I used:
copd.id <- DATA[!duplicated(DATA$patnr),]
copd.id$event <- as.numeric(!is.na(copd.id$deathdate))
lmeFit.copd <- lme(CysC~obstime+obstime:sex, random=~obstime|patnr, data=DATA)
coxFit.copd <- coxph(Surv(time,event)~sex, data=copd.id, x=TRUE)
jointFit.copd <- jointModel(lmeFit.copd, coxFit.copd, timeVar="obstime",method="piecewise-PH-aGH")
summary(jointFit.copd)
And I get the following error message:
jointFit.copd <- jointModel(lmeFit.copd, coxFit.copd, timeVar="obstime",method="piecewise-PH-aGH") Fehler in jointModel(lmeFit.copd, coxFit.copd, timeVar = "obstime", method = "piecewise-PH-aGH") : it seems that there are longitudinal measurements taken after the event times for some subjects (i.e., check subject(s): 0010000343, 0010000695, 0010000873, 0010002412, 0010002782, 0010003305, 0010003865, 0010003975, 0010004179, 0010004534, 0010004943, 0010005724, 0010007075, 0010007495, 0010008083, 0010008279, 0010008488, 0010008692, 0010008751, 0010009439, 0010010330, 0010011663, 0010012262, 0010012543, 0010012575, 0010013477, 0010014195, 0010015876, 0010016684, 0010017677, 0010018443, 0010019213, 0010019403, 0010019646, 0010020446, 0010020695, 0010021115, 0010021159, 0010021916, 0010022698, 0010024937, 0010026652, 0010030656, 0010031115, 0010031654, 0010031760, 0010033685, 0010034046, 0010034303, 0010035140, 0010037655, 0010038043, 0010038117, 0010038168, 0010038622, 0010039907, 0010042346, 0010044178, 0010046528, 0010046756, 0010048385, 0010049308, 0010049625, 0010049854, 0010050309, 0010051869, 0010052193, 0010052645, 0010052927, 0010053024, 0010054182, 0010055882, 001
summary(jointFit.copd) Fehler in summary(jointFit.copd) : Objekt 'jointFit.copd' nicht gefunden
The thing is: I checked the data, and there are no measurements after the events.
What am I missing here?
Try sorting DATA on patnr before running the models. I can induce the error you get for your DATA in the aids example by making a newid that is not sorted.
library(JM)
is.unsorted(aids$patient)
length(unique(aids$patient))
## make newid that is unsorted
no.rows <- rle(c(aids$patient))
aids$newid <- paste0(rep( rep(LETTERS, length.out=467), no.rows$lengths), aids$patient)
is.unsorted(aids$newid)
length(unique(aids$newid))
## change `patient` to newid in the following lines to get
## error message for `jointFit.aids`
aids.id <- aids[!duplicated(aids$patient),]
lmeFit.aids <- lme(CD4~obstime + obstime:drug, random=~obstime|newid, data=aids)
coxFit.aids <- coxph(Surv(Time,death)~drug,data=aids.id, x=TRUE)
jointFit.aids <- jointModel(lmeFit.aids, coxFit.aids, timeVar="obstime",method="piecewise-PH-aGH")
summary(jointFit.aids)
This came to mind when I looked at the code for JM on github. Since the error message is printed at line 139 in jointModel.R, I looked at how things were defined above that line and came to the conclusion that assumptions are made about the dataset being sorted on the id variable. So while it is true that no subject in DATA violates the condition, the order gets scrambled up in the code. In the absence of a reproducible example, my first hint was that the first 20 lines of the dataset given in the question had IDs 0010000343, 0010000873, 0010002412, but the error message posted lists the ids in sorted order had 0010000343, 0010000695, 0010000873, 0010002412 -- indicating that DATA is not sorted on the id variable (patnr) that is listed in the lme() model. The original aids example is sorted on id (patient).