drizopoulos / jm Goto Github PK
View Code? Open in Web Editor NEWJoint Models for Longitudinal & Survival Data under Maximum Likelihood
Joint Models for Longitudinal & Survival Data under Maximum Likelihood
Hi Dimitris,
I wondered why you did the change to the 'log.posterior.b' function. I met issues when running survfitJM and found that this was due to the change of 'log.posterior.b' function. Could you check this?
Many thanks,
Simon
Dear prof. Rizopoulos,
My spline-based JM (see below) runs but gives the warning mentioned in the title, and produces incorrect values, and NaN standard-errors, z-values and p-values.
From forums I read a solution that I do not understand?
https://stackoverflow.com/questions/53130611/error-in-nonlinear-optimization-problem-infinite-or-missing-values-in-x
surv_splines <- coxph(Surv(years, binstatus) ~ cluster(id) + disease+ age + gender + country, data = survtrain, x = TRUE, model=TRUE)
spline_lme <- lme(score~ ns(years, df=3) * (disease+ age + gender + country), na.action = na.exclude, random= list(id = pdDiag(form= ~ns(years, df=3))), data=longtrain
JMsplines <- jointModel(lme_splines_df3,spline_lme, timeVar = "years", method = "spline-PH-aGH")
I'm running the following in a dataset with 2641 subjects with follow-up in years and time to event between 0.01 and 4.9 years. The lme and the coxfit run with no problems but the jointModel returns
Error in aeqSurv(Y) : aeqSurv exception, an interval has effective length 0
Time-to-event
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.01369 1.14716 2.00685 2.05442 2.71869 4.86790
lme.NACC<-lme(Y ~ yrs + AGE + SEX + EDUC,
random = ~ yrs| NACCID, data=longdata_complete3)
coxfit.NACC<-coxph(Surv(time=timetoevent, event=event) ~ AGE + SEX + EDUC,
data=survdata_complete2, x=TRUE)
jointfit.NACC<-jointModel(lme.NACC,coxfit.NACC, timeVar="yrs",
method = "piecewise-PH-aGH")
Hi Dr. Rizopoulos,
I am having trouble getting the standard normal density function from your book to work with the DerivSplines and was wondering if I need to adapt it in some way. I am currently using
<
g <- function(u, pow = 0) {
f <- function(t) integrate(function(s) s^pow * dnorm(t - s), 0, t)$value
sapply(u, f)
}
iform_cew<-list(fixed=-1+I(g(time))+-1+I(g(time))+I(g(time,1)),indRandom=1:2)
ins(time,knots=3,weight.fun =g)+
I(g(time)*(cov)),
indFixed=1:5,
random=
corresponding to lme object with code
<
lmefit<-lme(longvar~ns(time,3)+cov,random=~time|ID,
data=data,control=lmeControl(opt="optim"))
Running the model with this derivForm returns an error: "Error in integrate(function(s) s^pow * dnorm(t - s), 0, t) :
evaluation of function gave a result of wrong length
In addition: Warning messages:
1: In s^pow :
longer object length is not a multiple of shorter object length
2: In s^pow * dnorm(t - s) :
longer object length is not a multiple of shorter object length"
I haven't been able to find any sample syntax for applying weighting functions to splines for JM and would appreciate any guidance. Thank you!
... I thought. But it does not. Everything is neat and fine with jointModel().
issue can be closed.
The issue was that lme() stores the name of the formula passed to it. Not the formula object itself. I think this might be a bug in lme().
When passing a list of formulas programmatically to lme() then the reference between the "call" and the "data" objects is broken.
jointModel() references correctly between call and data, but since this does not make any sense, jointModel() throws and error...
Hello there. Just starting with JM and having an issue that seems to be related to the size of the matrixes from the two sbumodels, even though they originate from the same df.
Your help would be greatly appreciated.
Here's the error that I get after running JMbayes2:
Error in jointModelBayes(mod_fi_onegroup, cox_fit, timeVar = "week") :
sample sizes in the longitudinal and event processes differ.
this is the reprex() {though it differs from the error message above}
#longitudinal submodel
mod_fi_onegroup=lme(fi ~ week + one_group,
random= ~ week|id,
corr=corAR1(form = ~ 1|id),
control=lmeControl(opt='optim'),
data=fi.long, method="REML",
na.action = na.exclude)
#> Error in lme(fi ~ week + one_group, random = ~week | id, corr = corAR1(form = ~1 | : could not find function "lme"
#coxph time-to-event submodel
cox_fit<-coxph(Surv(age_weeks, censored) ~ one_group + cluster (id),
data = fi.long, x=TRUE, model = TRUE)
#> Error in coxph(Surv(age_weeks, censored) ~ one_group + cluster(id), data = fi.long, : could not find function "coxph"
library(JMbayes2)
#> Warning: package 'JMbayes2' was built under R version 4.3.2
#> Loading required package: survival
#> Loading required package: nlme
#> Loading required package: GLMMadaptive
#> Warning: package 'GLMMadaptive' was built under R version 4.3.2
#> Loading required package: splines
joint_ca <- jm(cox_fit, mod_fi_onegroup, timeVar = "week")
#> Warning in jm(cox_fit, mod_fi_onegroup, timeVar = "week"): unknown names in
#> control: timeVar
#> Error in eval(expr, envir, enclos): object 'mod_fi_onegroup' not found
Originally posted by @marazzoli in #38 (comment)
As said, survfitJM
expects the data frame newdata
to be ordered by subject, and by time within each subject, else strange things happen.
This is not documented, as far as I could find.
An example based on chapter 7 of Joint Models for Longitudinal and Time-to-Event Data with Applications in R:
The code:
# Load packages JM and lattice
library("JM")
library("lattice")
set.seed(123) # we set the seed for reproducibility, early to get reproducible scrambling
# scramble rows
pbc2 <- pbc2[sample(seq_len(NROW(pbc2))), ]
# indicator for the composite event for the PBC dataset
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
# indicator for abnormal prothrombin time
pbc2.id$Pro <- with(pbc2.id, factor(prothrombin >= 10 & prothrombin <= 13,
labels = c("Abnormal", "Normal")))
pbc2$Pro <- rep(pbc2.id$Pro, tapply(pbc2$id, pbc2$id, length))
#################
# Section 7.1.3 #
#################
lmeFitBsp.pbc <- lme(
fixed = log(serBilir) ~ bs(year, 4, Boundary.knots = c(0, 15)),
random = list(
id = pdDiag(form = ~ bs(year, 4, Boundary.knots = c(0, 15)))),
data = pbc2)
coxFit.pbc <- coxph(Surv(years, status2) ~ drug + Pro,
data = pbc2.id, x = TRUE)
jointFitBsp.pbc <- jointModel(lmeFitBsp.pbc, coxFit.pbc,
timeVar = "year", method = "piecewise-PH-aGH")
# conditional survival probabilities for Patient 2 using Monte Carlo
survPrbs <- survfitJM(jointFitBsp.pbc, newdata = pbc2[pbc2$id == 2, ])
survPrbs
survPrbs$last.time[[1]] # 0
max(survPrbs$obs.times[[1]]) # but actually 8.8
plot(survPrbs, include.y = TRUE) # See image below, longitudinal outcome goes through survival plot
# estimate for longitudinal outcome is a line that goes back and forth through time, crossing itself
# When not including y, it may not be obvious that these predictions are incorrect
And when we estimate for multiple subjects, it goes even more wrong
survPrbsAll <- survfitJM(jointFitBsp.pbc, newdata = pbc2[pbc2$id %in% 1:5, ])
length(survPrbsAll$fitted.y[[1]]) # 7 ys fitted for first subject
length(survPrbsAll$fitted.times[[1]]) # However, these are fitted on 3 timepoints
plot(survPrbs, include.y = TRUE) # And we can't plot that
A workaround would be to set the order of newdata
within survfitJM.jointModel
:
newdata <- newdata[order(newdata[[idVar]], newdata[[object$timeVar]]), ]
Or to check order and error if it's not ordered properly:
stopifnot("newdata must be ordered by subject, and by time within each subject" = identical(order(newdata[[idVar]], newdata[[object$timeVar]]), seq_len(NROW(newdata)))
Alternatively, the code could be rewritten to no longer depend on ordering, but that's out of my league.
I can submit a single-line pull request for either of these fixes if that helps, however, I'm unsure if other functions are affected by the same behaviour.
I try to run sample codes from http://jmr.r-forge.r-project.org/Chapter5.html
pbc2.idCR <- crLong(pbc2.id, statusVar = "status",
censLevel = "alive", nameStrata = "CR")
lmeFit.pbc <- lme(log(serBilir) ~ drug * (year + I(year^2)),
random = ~ year + I(year^2) | id, data = pbc2)
coxFit4.pbc <- coxph(Surv(years, status2) ~ (drug + age) * CR + strata(CR),data = pbc2.idCR, x = TRUE)
jointFit11.pbc <- jointModel(lmeFit.pbc, coxFit4.pbc,
timeVar = "year", method = "spline-PH-aGH", CompRisk = TRUE,
interFact = list(value = ~ CR, data = pbc2.idCR))
lme and coxph ran fine but jointModel produced an error as following
'Error in exp(coefs[1]) : non-numeric argument to mathematical function'
environment:
macOS Big Sur (version 11.3.1)
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Platform: x86_64-apple-darwin17.0 (64-bit)
Could you please help
Thanks!
Hi Prof Rizopoulos,
Is there an implementation that would allow an entirely new subject (outside the model training data) to be used to create dynamic predictions? If for example we have their baseline covariates and longitudinal biomarker data are known.
Thank you in advance.
Hi @drizopoulos,
I am fitting some joint models of organ dysfunction in critical illness. In order to help with stability, I center and scale variables before modelling (including the outcome of the longitudinal submodel). Some of the values of the longitudinal outcome will now take negative values. There is a strange appearance to the ROC plots (see reproducible example from the package sample data) where the curves don't complete their path. The AUC results are also much lower than would be expected. I can't seem to isolate the issue to any other cause. Is this normal behaviour? Must variables inside the joint model only take positive values? is it necessary that other variables in the model are also positive (e.g. explanatory variables inside the longitudinal or cox submodels)?
library(JM)
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
# center and scale longitudinal outcome
bil_scaling <- scale(pbc2$serBilir)
# reassign to data frame
pbc2$serBilir <- bil_scaling[, 1]
pbc2.id$serBilir <- pbc2$serBilir[pbc2$year == 0]
# fit models
lmeFit <- lme(serBilir ~ ns(year, 3),
random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "year",
method = "piecewise-PH-aGH")
# ROC plots
ND <- pbc2[pbc2$id == "7", ]
roc1 <- rocJM(jointFit, dt = c(2, 4, 8), data = ND, idVar = "id")
plot(roc1)
I couldn't find any other reference to this in your book or on stack overflow. Thank you in advance for your time and insights on this.
Dear Dimitris,
I'm trying to use the JM library, but I get error messages and they seem connected to the time duration units: days or months.
I was using months as time scale and I then got the message:
"Error in aeqSurv(Y) :
aeqSurv exception, an interval has effective length 0"
It seems that this happens if I use months for survival analysis.
If I change my time variables into duration in days, I get another error message:
"Error in if (t1 || t2) { : missing value where TRUE/FALSE needed"
I found a similar issue here in the GitHub and the advice was not to use days as unit for duration.
I'm a bit stuck with this as nor the days neither the months seem to work. I would very much appreciate it if you could give me some advice?
The code for Figure 5-1 in the book "Joint Models for Longitudinal and Survival Data" on page 103 is as follows
Dear,
Thank you again for this software. Unfortunately I am having trouble getting it to run on my own data. When I run the code below , I get the following error message from JointModel(): Error in solve.default(VC) : system is computationally singular: reciprocal condition number = 3.40196e-17 , Error: cannot allocate vector of size 86.8 Gb,
Error in lme.formula(obs ~ IND_SEVVP0 + AUDI0_C + SEXE + CENTRE + DIPNIV0C + :
nlminb issue, code d'convergence error = 1 message = false convergence (8)
Error in optim(thetas, opt.survWB, gr.survWB, method = "BFGS", control = list(maxit = if (it < : unfinished value provided by optim
Warning message:
In jointModel(lme_Isa_30_M4VP2, Cox_Isa_30_M4VPquad, timeVar = "T", :
infinite or missing values in Hessian at convergence.
MixteISAVP4 <- Base_ISA[(!is.na(Base_ISA$IND_SEVVP0) & (!is.na(Base_ISA$AUDI0_C)) & (!is.na(Base_ISA$DIPNIV0C))
& (!is.na(Base_ISA$VIVRE_SEUL0)) & (!is.na(Base_ISA$REVENU0C)) & (!is.na(Base_ISA$FUME0))
& (!is.na(Base_ISA$BMI0C)) & (!is.na(Base_ISA$ATCDCAR)) & (!is.na(Base_ISA$ATCDAVC))
& (!is.na(Base_ISA$HTA0_1)) & (!is.na(Base_ISA$DEPRES0C)) & (!is.na(Base_ISA$DIABBIS0C))
& (!is.na(Base_ISA$TRIGLY0C)) & (!is.na(Base_ISA$APOE4C)) & (!is.na(Base_ISA$HYPCT024C))), ]
CoxISA4$DELAIS<- CoxISA4$DELAIS/365
MixteISAVP4$T <- MixteISAVP4$T/365
ctrl <- lmeControl(opt='optim');
lme_Isa_30_M4VP2 <- lme(obs ~ IND_SEVVP0 + AUDI0_C + SEXE + CENTRE + DIPNIV0C + Prem + Age65 + VIVRE_SEUL0 + REVENU0C +
FUME0 + BMI0C + ATCDAVC + ATCDCAR + HTA0_1 + DEPRES0C + DIABBIS0C + TRIGLY0C +
APOE4C + HYPCT024C + T + I(T^2) +
IND_SEVVP0T + AUDI0_CT + SEXET + CENTRET + DIPNIV0CT + Age65T + DIABBIS0CT + APOE4CT +
IND_SEVVP0I(T^2) + AUDI0_CI(T^2) + SEXEI(T^2) + CENTREI(T^2) + DIPNIV0CI(T^2) + Age65I(T^2) + DIABBIS0CI(T^2) + APOE4CI(T^2),
random = ~ T + I(T^2) | NUM, method="ML", control=ctrl, na.action=na.omit, data = MixteISAVP4)
#summary(lme_Isa_30_M4VP2)
CoxISA4 <- MixteISAVP4[!duplicated(MixteISAVP4$NUM), ] # passe en 1 ligne par sujet
Cox_Isa_30_M4VPquad <- coxph(Surv(DELAIS, INDICDP) ~ IND_SEVVP0 + AUDI0_C + SEXE + CENTRE + DIPNIV0C + Prem + Age65
+ VIVRE_SEUL0 + REVENU0C + FUME0 + BMI0C + ATCDAVC + ATCDCAR + HTA0_1 + DEPRES0C
+ DIABBIS0C + TRIGLY0C + APOE4C + HYPCT024C, data = CoxISA4, x = TRUE, model=true)
#Cox_Isa_30_M4VPquad
ctrljm<-list(iter.EM=500)
fitJOINTvalue_ISA30_M4VPquad <- jointModel(lme_Isa_30_M4VP2,
Cox_Isa_30_M4VPquad,
timeVar = "T",
method="Cox-PH-GH",
verbose=TRUE,
control=ctrljm)
#control = list(GHk = 3, lng.in.kn=1))
fitJOINTvalue_ISA30_M4VPquad <- jointModel(lme_Isa_30_M4VP2,
Cox_Isa_30_M4VPquad,
timeVar = "T",
method = "spline-PH-aGH")
fitJOINTvalue_ISA30_M4VPquad <- jointModel(lme_Isa_30_M4VP2,
Cox_Isa_30_M4VPquad,
timeVar = "T",
method = "piecewise-PH-aGH")
summary(fitJOINTvalue_ISA30_M4VPquad)
jointFit.aids3 <- jointModel(lme_Isa_30_M4VP2, Cox_Isa_30_M4VPquad,
timeVar = "T", method = "piecewise-PH-aGH", GHk = 3)
jointFit.aids6 <- jointModel(lme_Isa_30_M4VP2, Cox_Isa_30_M4VPquad,
timeVar = "T", method = "piecewise-PH-aGH", GHk = 6)
jointFit.aids9 <- jointModel(lme_Isa_30_M4VP2, Cox_Isa_30_M4VPquad,
timeVar = "T", method = "piecewise-PH-aGH", GHk = 9)
Hi, @drizopoulos,
I am trying to use derivForm with splines and getting the following error:
"Error in Xderiv %*% fixef(lmeObject)[derivForm$indFixed] :
non-conformable arguments"
Can anyone show an example how to set the derivForm when you use splines in the linear mixed effect model?
Here is my code
fitLME_poly <- lme(sofa_recent ~ ns(hr_day, 3)*ethnicity+U_sofa_adm+age10+male,
random = list(stay_id = pdDiag(form = ~ ns(hr_day, 3))), data = mdat)
fitCox_poly <- coxph(Surv(icu_day,died) ~ ethnicity+age10+male+U_sofa_adm,data=mdat.id,x=TRUE)
dForm <- list(fixed = ~ dns(hr_day, 3) + dns(hr_day, 3)*ethnicity, indFixed = c(2,3,4,11,12,13,14,15,16,17,18,19),
random = ~ dns(hr_day, 3) , indRandom = c(2,3,4))
fitJM_pwph_poly <- jointModel(fitLME_poly,fitCox_poly,timeVar="hr_day",
method="piecewise-PH-aGH",
parameterization = "both", derivForm = dForm,
control=list(verbose=TRUE,iter.EM=100))
Thanks
I have encountered the caught segfault error,
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
*** caught segfault ***
address (nil), cause 'unknown'
*** recursive gc invocation
*** recursive gc invocation
and even with tryCatch
, the error cannot be skipped (It seems that the segfault error is not as simple as the errors can be handled with try exception?)
I can reproduce it with the following simulation setting,
set.seed(1234)
n = 100
time2event = sample(1:3, n, replace = TRUE) + 100
censor = numeric(n)
censor[time2event == 103] = 1
treatment = sample(1:4, n, replace = TRUE)
covariate = rnorm(n)
t = runif(n)
id = 1:n
df = data.frame(id, time2event, censor, treatment, t, covariate)
cox.fit = coxph(Surv(time2event, censor) ~ treatment, data = df, x = TRUE)
the above setting is to construct a bad survival model, where Concordance= NaN (se = NaN )
,
> summary(cox.fit)
Call:
coxph(formula = Surv(time2event, censor) ~ treatment, data = df,
x = TRUE)
n= 100, number of events= 31
coef exp(coef) se(coef) z Pr(>|z|)
treatment 3.061e-17 1.000e+00 1.905e-01 0 1
exp(coef) exp(-coef) lower .95 upper .95
treatment 1 1 0.6884 1.453
Concordance= NaN (se = NaN )
Likelihood ratio test= 0 on 1 df, p=1
Wald test = 0 on 1 df, p=1
Score (logrank) test = 0 on 1 df, p=1
then run
lme.fit = lme(covariate ~ t, random = ~ 1| id, data=df)
for (i in 1:100) {
try({
print(i)
jointModel(lme.fit, cox.fit, timeVar = "t")
})
}
it throws
[1] 1
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
[1] 2
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
[1] 3
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
...
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
[1] 26
Error in names(fit$coef) <- cname :
'names' attribute [4] must be the same length as the vector [2]
[1] 27
*** caught segfault ***
address (nil), cause 'unknown'
*** recursive gc invocation
*** recursive gc invocation
note that the segfault
is not immediately thrown, and instead at iter=27, and once it is thrown, the for loop also exit, even though I have wrapped the command with try
.
I am using JM_1.5-2 in R4.2.0, and here is the detailed session info,
> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.6 (Maipo)
Matrix products: default
BLAS/LAPACK: /home/*****/miniconda3/lib/libopenblasp-r0.3.21.so
locale:
[1] C
attached base packages:
[1] splines stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] JM_1.5-2 survival_3.4-0 nlme_3.1-160 MASS_7.3-58.1
loaded via a namespace (and not attached):
[1] compiler_4.2.0 Matrix_1.5-1 grid_4.2.0 lattice_0.20-45
Hi Dr. Rizopoulos,
Like a few of the commentors below, I am having an issue with getting my joint model to run, receiving an error
"Error in aeqSurv(Y) :
aeqSurv exception, an interval has effective length 0"
fit.lm1<-lme(TSH_~ time1,random= ~time1 | code,data=dataLong)
model.cox<-coxph(Surv(time_event_yr,DM4)~BMI_3,data=data,x=T,model=TRUE)
summary(model.cox)
aeqSurv(model.cox, tolerance = sqrt(.Machine$double.eps)
Error in aeqSurv(model.cox, tolerance = sqrt(.Machine$double.eps)) :
argument is not a Surv object
Any advice for how to get my model to run would be very appreciated!
Thanks!
Hi everyone,
I have a couple of questions about what the simulate function is capable of.
Can it be utilized to sample out the evolution (paid and time) of open claims (unevolved data points) using N simulations?
In addition, how do the glmmPQL and LME functions interact between each other, right now they are probably equivalent (or not?) as I am using gaussian family, but would like to try to work using other dist families (gamma, inv gauss) - currently these don't seem to be working.
Thanks
John Gribble
Can you extract or create the lme object from the jointModel function after you have fit the JM? If so, how do you go about doing it?
Hi Dr. Rizopoulos,
Like a few of the commentors below, I am having an issue with getting my joint model to run, receiving an error
fit.lm1<-lme(log(TSH_)~time ,random= ~time | code,control=lmeControl(opt="optim"),data=dataLong)
summary(fit.lm1)
model.cox<-coxph(Surv(Time,DM1)~sex15_1,data=data,x=T,model=TRUE)
summary(model.cox)
fit.JM<-jointModel(fit.lm1,model.cox,timeVar="time",method ="piecewise-PH-aGH")
Error in d * log.hazard + log.survival : non-conformable arrays
In addition: Warning message:
In wk * rep(x$P, each = nk) :
longer object length is not a multiple of shorter object length
Any advice for how to get my model to run would be very appreciated!
Thanks
Hi sorry if I am missing something obvious but I get an error when I try and use the basic diagnostic plot for JM if I just have a random intercepts model. So in terms of an example you used in "JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data".
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ 1 | patient, data = aids)
fitSURV <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
fit.JM <- jointModel(fitLME, fitSURV, timeVar = "obstime", method = "piecewise-PH-aGH")
plot(fit.JM, which = 1)
plot(fit.JM, which = 2)
plot(fit.JM, which = 3)
plot(fit.JM, which = 4)
Here I have just changed the random effects in the mixed effects model to 1|patient instead of obstime|patient and it comes up with the following error when I try and plot the 3rd and 4th plot.
Error in Zs * b[id.GK, , drop = FALSE] : non-conformable arrays
If you could help me understand why that is it would be greatly appreciated! again sorry if this is obvious :)
I don't currently have a reproducible example in hand, though I'll work on one.
I'm running into this error message having called jointModel
,
Error in aeqSurv(Y) :
aeqSurv exception, an interval has effective length 0
I tracked the aeqSurv
call down to coxph
here. I noticed really large values in the Time
variable passed to initial.surv
. I tracked these down to this spot.
I'm wondering about the reason for the exp
here?
Time <- exp(survObject$y[, 1])
Dear Professor Rizopoulos,
I'm running through the section on "parameterization of random effects". And I'm confused about how to calculate the standard deviation when correcting for a particular subject's slope by dividing the original value by its corresponding standard deviation.
I appreciate your reply!
I got an error about some observations in my data having measurements on different temporal ranges. After some investigation I found that this was due to a comparison of floating point numbers. I have submitted a pull request that fixes it (2e5ec8d).
Are you considering expanding the capabilities of JM to accept mixed models fitted by the lme4
package, in addition to those from nlme
? It seems that lme4
is still under active development and it's often faster, too, especially with 'crossed' random effects. Thanks for a very useful package, by the way!
I am getting an error message when attempting to extract event model residuals for competing risks joint models (JM version 1.4-5). A reproducible example is below.
Fit the competing risks model from the jointModel
help documentation
library(JM)
data("pbc2.id")
data("pbc2")
# Fit model from jointModel help documentation
pbc2.idCR <- crLong(pbc2.id, "status", "alive")
lmeFit.pbc <- lme(log(serBilir) ~ drug * year,
random = ~ year | id,
data = pbc2)
coxCRFit.pbc <- coxph(Surv(years, status2) ~ (drug + sex)*strata + strata(strata),
data = pbc2.idCR,
x = TRUE)
jmCRFit.pbc <- jointModel(lmeFit.pbc, coxCRFit.pbc, timeVar = "year",
method = "spline-PH-aGH",
interFact = list(value = ~ strata, data = pbc2.idCR),
CompRisk = TRUE)
Then, if I run
residuals(jmCRFit.pbc, process = "Event")
the following error message is received: Error in rep(object$y$d, ni) : invalid 'times' argument
.
Interestingly, if I run
residuals(jmCRFit.pbc, type = "Martingale")
i.e. without any process
argument, I get output. However, I am not sure what the output is.
If these residuals are not supposed to be accessed using competing risks models, perhaps adding a if (compRisk) stop() else ...
line into the code would be helpful, along the same lines as was done for when a user tries to use MI = TRUE
with such models.
Dear prof. Rizopoulos,
I want to calculate survival probabilities in a 100-patient sample.
All JM functions work lovely, except surfitJM, which throws the abovementioned error.
What does it imply?
survfitJM <- JM:::survfitJM.jointModel
sprobs <- survfitJM(JM,testdata,idVar="id",survTimes=c(28/365.25,90/365.25))
Error in rowMeans(rr, na.rm = TRUE) : 'x' must be numeric
I have a joint model with the form
jointFit<- jointModel(lmeFitBsp.pro, coxFit.pro,
timeVar = "time", method = "piecewise-PH-aGH",
knots=c(150,175,200,225,250,260,270,280,300))
where the events come much later in the time period so the default 7 knots will not work and I am setting them manually
when I try to get a roc :
rocJM(jointFit, dt = 5, data = plcbData,
+ M = 1000, burn.in = 500)
I get
Error in optim(rep(0, ncz), ff, gg, y = y.new, tt = survMats.last, mm = method, :
initial value in 'vmmin' is not finite
In addition: There were 50 or more warnings (use warnings() to see the first 50)
Warning messages:
In wk * rep(P, each = 7) :
longer object length is not a multiple of shorter object length
I think it is realted to this where the n = 7 default knots is hard coded:
Line 62 in 05bfde1
Hi Dr. Rizopoulos,
In "survfitJM.jointModel.R" line 180, you may forget to add "Dalpha.new <- thetas.new$Dalpha".
BTW, I have another question about line 250 in "splinePHGH.fit.R". Based on the closed formula of sigma.hat, it should be the trace of Z'Zvb, why directly summing the matrix instead of summing the diagonal of the matrix?
Thanks,
Xiangyu
Hi everyone,
I had this error when running fuction jointModel() in JM package with method = "Cox-PH-aGH" option.
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
pbc2$lserBilir <- log(pbc2$serBilir)
lmeFit1 <- lme(lserBilir ~ year + drug:year, data = pbc2,
random = ~ year | id)
survFit1 <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit1 <- jointModel(lmeFit1, survFit1, timeVar = "year",
method = "Cox-PH-aGH")
Error in model.frame.default(data = FALSE, formula = ~long + W1) :
'data' must be a data.frame, environment, or list
I cannot fix it since pbc2 and pbc2.id are originally data.frame object.
It runs except with this method option.
Does anyone have any clue?
Thanks!
from the book:
pbc <- pbc2[c("id", "serBilir", "drug", "year", "years",
"status2", "spiders")]
pbc$start <- pbc$year
splitID <- split(pbc[c("start", "years")], pbc$id)
pbc$stop <- unlist(lapply(splitID,
function (d) c(d$start[-1], d$years[1]) ))
pbc$event <- with(pbc, ave(status2, id,
FUN = function (x) c(rep(0, length(x)-1), x[1]) ))
tdCox.pbc <- coxph(Surv(start, stop, event) ~ drug * spiders + cluster(id),
data = pbc, x = TRUE, model = TRUE)
jointFit8.pbc <- jointModel(lmeFit.pbc, tdCox.pbc, timeVar = "year",
method = "spline-PH-aGH")
Error in jointModel(lmeFit.pbc, tdCox.pbc, timeVar = "year", method = "spline-PH-aGH") :
use argument 'model = TRUE' and cluster() in coxph().
the code test for the cluster slot but it is now (cluster)
Hi Dr. Rizopoulos,
I am slightly concerned with the point estimates for my association parameters for a current value + current slopes joint model. I am using a large dataset with 91000 observations from 3650 participants. Within a piecewise-PH-GH formulation, estimates for my association parameters for current values plus slope seem funky compared to other parameterizations, and across different baseline hazards (piecewise vs. weibull) they are completely different.
When I used method= "piecewise-PH-GH" in the current values and cumulative effects parameterizations, I got estimates of assoct = 0.2 (SE = 0.01) and assoct.s= 0.01 (SE =0.001) , respectively. When using the current values+ current slopes parameterizations I got estimates of assoct. = -2.05 (SE= 0.25) and assoct.s = 2.25 (SE= 0.25). Would this sort of difference in hazard ratios be plausible in absence of some sort of issue in numerical integration?
When I run the current values+ current slopes using method= "weibull-PH-GH", I get point estimates of assoct = 0.36 (SE = 0.24) and assoct.s = -0.16 (SE=0.24). Is it typical for point estimates to change like this when using different baseline hazard formulations? Is there a way to tell which one is more correct?
Thanks for all the help!
I got a warning after running predict with a jointModel
for a second time:
Warning message:
In rm(list = ".Random.seed", envir = globalenv()) :
object '.Random.seed' not found
The reason was that there was no seed to delete with rm
as it was deleted in the last run and no new seed was created.
I have submitted a pull request (#6) in which I moved the line
rm(list = ".Random.seed", envir = globalenv())
inside the else
clause.
(For a moment I thought it was an error and created the PR, otherwise I may not have bothered you as this does not affect the results.)
Hi, thanks for the amazing package. I am wondering in the two parts of the joint modeling framework, can they accept patients entering the study/observation from different time points?
I am dealing with data from a patient database, therefore, in a real-world setting, patients come to visit the doctor and with their blood sample drawn at totally different time points/intervals. Some patients started early some with their initial recording later.
So my question is whether the survival part of the Joint Modeling within this package can be fitted using the following form?
coxFit.DF <- coxph(Surv(start, stop, event) ~ sex + age + group, data, x = TRUE)
Because when I fitted the above Cox model with the difference between the times they were followed time = stop-start
:
coxFit.DF <- coxph(Surv(time, event) ~ sex + age + group, data, x = TRUE)
There were errors telling me that some patients had observations after the events. And when I use the Surv(start, stop, event)
one, error message also came up and ask me to refit the Cox model using x = TRUE
argument which I had already added.
One of the flexibility of the time-dependent Cox regression model is that patients are allowed to enter the study from different start points, If there is also a way to allow for this in the Joint Modeling framework, I would like to see how and it will be greatly appreciated.
Many thanks.
Dear prof,
In a large dataset (8081 patients, 69594 measurements) aucJM (JM:::aucJM.jointModel) keeps giving the error below, but if I change Tstart or Tstop slightly, it does work. My data is complete, sorted on id and time of measurement, I have no events on t=0.
Why does it give this error?
AUC90d_base <- aucJM(JM_spline_valslope, longtestsurv, Tstart = 0.0001, Thoriz = 92/365.25, idVar = "id")
Error in FUN(X[[i]], ...) : subscript out of bounds
traceback()
3: lapply(X = X, FUN = FUN, ...)
2: sapply(pi2$summaries, "[", 1, 2)
1: aucJM(JM_spline_valslope, longtestsurv, Tstart = 1e-04, Thoriz = 92/365.25,
idVar = "id")
The JM is speficied as:
survFit <- coxph(Surv(tyrs, binstatus) ~ age10 + femalegender + cluster(id),
data = survtrain, x = TRUE)
lme_splines <- lme(score~ ns(labyrs, df=3)* (age10 + femalegender),
random= list(id = pdDiag(form= ~ns(labyrs, df=3))), data=longtrain)
dform = list(fixed = ~ 0 + dns(labyrs, df=3) + dns(labyrs,df=3):age10 +dns(labyrs,df=3):femalegender,
random = ~ 0 + dns(labyrs, df=3),
indFixed = c(2:4,7:12), indRandom = c(2:4))
#JM
JM_spline_valslope <- jointModel(lme_splines, survFit, timeVar = "labyrs",
method = "spline-PH-aGH", param = "both", derivForm = dform)
Hello Dr Rizopoulos,
I am simulating data coming from an Exponential PH model and I am trying to feed it to the jointModel() function to see whether I get back my parameter estimates.
You have mentioned in page 95 of your joint modelling book that survreg fits AFT model and thus if we want MLE's under the relative risk parameterization we will have to supply relevant initiation of the values to the jointModel() function.
In specifying the "weibull-PH-aGH" , I noticed that my association parameter was extremely far from what I had simulated (ie. -0.0006 vs 0.2). However, without using survreg (using coxph) and calling jointModel directly (without init), I was able to get back my association parameter (0.1957).
I am wondering what did I do wrong or if theres something I should be aware of?
fitLME<- lme(measurements ~ Time + FixedCovariate, random= ~ Time|ID, data = longitudinalDataframe)
marginal.survJM <- survreg(Surv(times, status)~W1+W2+X2, data=survivalDF, dist ="exponential", x=TRUE) #this is the Exponential AFT fit
### CONVERT PARAMETERS AS PER PAGE 95 of Dr Rizopoulos' book.
init.list<- list(betas = fixef(fitLME), sigma = fitLME$sigma, D = getVarCov(fitLME),
gammas = -coef(marginal.survJM)/marginal.survJM$scale,
sigma.t = 1)
JMFITWB<- jointModel(fitLME, marginal.survJM, timeVar= "Time", scaleWB=1, init = init.list, method="weibull-PH-aGH")
summary(JMFITWB)
coxPhModel <- coxph(Surv(times, status)~W1+W2+X2, data=survivalDF, x=TRUE)
JMFITCOX <- jointModel(fitLME, coxPhModel, timeVar = "Time", scaleWB=1)
summary(JMFITCOX)
Really would appreciate any help or guidance.
Thank you,
Faith
Hi,
I am a beginner in the area of Joint Models in survival analysis. While running R codes for Joint Model for Longitudinal and time to event data, I got the following error:
**Error in optim(thetas, opt.survPC, gr.survPC, method = "BFGS", control = list(maxit = if (it < :
non-finite value supplied by optim**
This is my code for linear mixed model, Survival model, and Joint model.
**> lmeFit<- lme (FUCEAvalue ~ year + SEX, random = ~ year| ID, data = longi_data.new)
coxFit<- coxph (Surv(year,DEATH)~ SEX, data = term.evnt.new, x="TRUE")
jointFit<- jointModel (lmeFit, coxFit, timeVar = "year", method = "piecewise-PH-aGH")**
Any information and advice you could give is much appreciated.
i am getting this error when i give the command for running a plspm
requesting someone to pls help me.
loyalty_pls1 = plspm(loyalty, loyalty_path, loyalty_blocks, modes = loyalty_modes)
Error in if (w_dif < specs$tol || iter == specs$maxiter) break :
missing value where TRUE/FALSE needed
I don't know why this error came up here. can somebody explain it to me.
loyalty_pls = plspm(loyalty, loyalty_path, loyalty_blocks, modes=loyalty_modes)
Error in if (w_dif < specs$tol || iter == specs$maxiter) break :
missing value where TRUE/FALSE needed
In addition: Warning message:
Dear Professor Rizopoulos,
I have some doubts regarding the interpretation of exp(Assoct) in joint model.
In your book "Joint models for longitudinal and time-to-event data _ with applications in R", section 4.1.1 The Survival Submodel mentions:
h0(t) exp{γ⊤wi + αmi(t)}
exp(α) denotes the relative increase in the risk for an event at time t that results from one unit increase in yi(t) at the same time point. `However, here mi(t) is defined as:
yi(t) = mi(t) + εi(t),
mi(t) = x ⊤ i (t)β + z ⊤ i (t)bi
the mi(t) was linear.
Because my y is non-linear, I used bs(time,3)
to fit the nonlinear process of y.
Below, I've included the code and results:
fitLME <- lme(y~ bs(chartday,3), random = ~ 1+chartday|subject_id,
data = data_for_lme)
fitSURV <- coxph(Surv(time, status) ~ apsiii+charlson_pasthistory_charlson_comorbidity_index,
data = data_for_cox, x = TRUE , model = TRUE)
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "chartday")
summary(fitJOINT)
Call:
jointModel(lmeObject = fitLME, survObject = fitSURV, timeVar = "chartday")
Data Descriptives:
Longitudinal Process Event Process
Number of Observations: 74816 Number of Events: 1391 (20.7%)
Number of Groups: 6713
Joint Model Summary:
Longitudinal Process: Linear mixed-effects model
Event Process: Weibull relative risk model
Parameterization: Time-dependent
log.Lik AIC BIC
-37348.72 74723.43 74811.99
Variance Components:
StdDev Corr
(Intercept) 0.6976 (Intr)
chartday 0.0807 -0.2029
Residual 0.2904
Coefficients:
Longitudinal Process
Value Std.Err z-value p-value
(Intercept) 5.0622 0.0088 572.5697 <0.0001
bs(chartday, 3)1 -0.6039 0.0163 -36.9510 <0.0001
bs(chartday, 3)2 1.5291 0.0348 43.9721 <0.0001
bs(chartday, 3)3 0.1008 0.0443 2.2755 0.0229
Event Process
Value Std.Err z-value p-value
(Intercept) -5.1218 0.1958 -26.1581 <0.0001
apsiii 0.0210 0.0009 24.0735 <0.0001
charlson_pasthistory_charlson_comorbidity_index 0.1153 0.0089 12.9207 <0.0001
Assoct -0.3082 0.0285 -10.8297 <0.0001
log(shape) 0.0302 0.0234 1.2906 0.1968
Scale: 1.0307
Integration:
method: (pseudo) adaptive Gauss-Hermite
quadrature points: 3
Optimization:
Convergence: 0
Can I still interpret exp(Assoct) as denoting the relative increase in the risk for an event at time t that results from one unit increase in yi(t) at the same time point at this point? Can exp(Assoct) be interpreted as this no matter how mi(t) is fitted (bs(), ns(), poly()...)?
Thank you very much!
Best,
Qian
Hi,
I am a beginner of Joint Models and JM package. While running jointModel() to fit joint models with competing risks, I got the following error:
Error in in optim(thetas, LogLik.splineGH, Score.splineGH, method = "BFGS", optim replies to the infinite value
My code as follow:
fitLME.null <- lme(score_L ~ followyear_L,
random = ~1 | id , data = data2, method = "REML")
coxCRFit.pbc<- coxph(Surv(totalfollow, status2) ~(age +BMI_L+hyp_L+exersice_L)*strata+strata(strata),data = data2.idCR,
x = TRUE)
jmCRFit.pbc <- jointModel(fitLME.null, coxCRFit.pbc, timeVar = "followyear_L",
method = "spline-PH-aGH",
interFact = list(value=~strata, data= data2.idCR),
CompRisk = T)
score_L is a scale score that ranges from 0 to 30. Status2 was the event including alive, disease and death.
All subjects had at least two follow-ups with id as personal identifier.
Any information and advice you could give is much appreciated.
I had this error when running
model.cox<-coxph(Surv(surv_time,incidenceDM)sex15_1+age+FHD+ waist+fbs+bs2HR+BMI,data=data.cox,x=T,model=TRUE) time,random= ~time | code,data=dataLong)
summary(model.cox)
fit.lm<-lme(logTSH
summary(fit.lm)
dform<-list(fixed=~1,inFixed=2,random=~1,indRandom=2)
fit.JM<-jointModel(fit.lm,model.cox,timeVar="time",parameterization="both",derivForm =dform)
summary(fit.JM)
Error in tapply (Long.deriv,id,tail,1):
Arguments must have same length
Thanks,
Hi Dr. Rizopoulos,
Like a few of the commentors below, I am having an issue with getting my joint model to run, receiving an error
"Error in aeqSurv(Y) :
aeqSurv exception, an interval has effective length 0"
My survival object uses the survreg function with weeks as my unit of time.
lmefit<-lme(marker~ns(week,3),random=~week|ID, data=Data,control=lmeControl(opt="optim")) survfit<-survreg(Surv(time,status)~covariate,data=Data.id,x=TRUE) fitjoint<-jointModel(lmefit,survfit,timeVar="week")
I haven't been able to get any of the fixes listed here to work. The issue doesn't seem to be one of 2 weeks vs 2.0 weeks. Entering "timefix=FALSE" into the survreg call gets an error for an unused argument, and adding
aeqSurv(survfit, tolerance = sqrt(.Machine$double.eps))
returns an error that survfit isn't a Surv object.
Any advice for how to get my model to run would be very appreciated!
Thanks!
Dear Dimitris Rizopoulos,
I try to use your JM package in a shiny app to plot the predicted conditional probabilities of survival with the plot of my longitudinal marker. However, my shiny app display an error when executing the following line of code :
plot.survfitJM(survPreds[[length(survPreds)]], estimator="median", conf.int=TRUE, include.y=TRUE)
survPreds being survfitJM object.
The displaying error is : comparison of these types is not implemented
Note that without the option include.y=TRUE it work fine.
To investigate the error I saved the R environment in the shiny app just before the problematic call to plot function. Then I opened the obtained .Rdata in a standard R and I manually execute the previous plot function. In thiw way I obtained the desired graph without error !
Alltogather, I have no idea where this error came from... could you help me ? If necessary I could provide you the entire code of my shiny app.
Thanks a lot.
Kind regards,
Florent
Hi,
I'm encountering an error while attempting joint modeling for my dataset using the spline-PH-GH method.
Here's a snippet of my code:
lmefit<-lme(y ~ time+age+sex+comorb+stage,random=~ time|ID,data = long)
coxfit<-coxph(Surv(time,DEATH)~age+sex+comorb+stage,data = Surv,x=T)
jointFit<-jointModel(lmefit, coxfit, timeVar = "time",method = "spline-PH-GH")
I got the following error message:
Error in if (t1 || t2) { : missing value where TRUE/FALSE needed.
This issue doesn't occur when I use other methods like "piecewise-PH-aGH" or "Weibull-AFT-GH".
I've attached a subset of my dataset for reference.
datasets.zip
I'd appreciate any insights into why this error is occurring and how I might resolve it.
Thanks in advance for your help!
Hi,
Thanks for the amazing work in JM. I wanted to fit a joint model with competing risks data using the specification of piecewise baseline hazard function. Does this package offer this option?
Thanks,
Jeffrey
Hi,
Thank you for this amazing software. Unfortunately I am having trouble getting it to run on my own data. When I run the code below (modelled on one of your examples) I get the following error message from JointModel():
Error in while (iter < maxits && !converged) { :
missing value where TRUE/FALSE needed
This is my code:
library(nlme)
library(survival)
library(JM)
data_dir <- "C:/mydir/"
cox_wide <- read.csv(paste(data_dir, "data_sm_wide.csv", sep=""), header=TRUE, sep=",")
cox_long <- read.csv(paste(data_dir, "data_sm.csv", sep=""), header=TRUE, sep=",")
ctrl <- lmeControl(opt='optim');
bloodfit <- lme(Age ~ tstart + tstart:Urate, random=~tstart|ID, method="ML", control=ctrl, data=cox_long, na.action=na.omit)
summary(bloodfit)
coxphobject <- coxph(Surv(time=duration, event=cstatus)~Urate, data=cox_wide, x=TRUE, model=TRUE)
summary(coxphobject)
# Joint Model
ctrljm<-list(iter.EM=500)
jmfit <- jointModel(bloodfit, coxphobject, method="Cox-PH-GH", timeVar = "tstart", verbose=T, control=ctrljm)
summary(jmfit)
I have attached a small subset of my data. It crashes after iteration 6. I am using R 3.5.1.
Thanks for your help.
Annette
data.zip
Hi, It looks like with the update that the random intercept diagnostic plots work but the random intercept and slope no longer works.
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ obstime | patient, data = aids)
fitSURV <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
fit.JM <- jointModel(fitLME, fitSURV, timeVar = "obstime", method = "piecewise-PH-aGH")
plot(fit.JM, which = 1)
plot(fit.JM, which = 2)
plot(fit.JM, which = 3)
plot(fit.JM, which = 4)
For the final two plots here I get the error
Error in Zs * b[id.GK, , drop = FALSE] : non-conformable arrays
Thank you for your time!
Hi Dr. Rizopoulos, thanks again for all the help!
I am trying to fit a cumulative effects model based on the following lme object
<lme_fit<-lme(var~ns(time,3)+cov,random=~time|ID,
data=data,control=lmeControl(opt="optim"))>
I haven't been able to get the derivForm argument right though. I've been trying
<iform<-list(fixed=-1+time+ins(time,3)+(cov*time),indFixed=c(1:5),-1+time+time^2/2,indRandom=c(1:2))>
random=
and getting an error "Error in terms.formula(derivForm$random) : invalid model formula in ExtractVars"
Would you be able to provide any guidance as to what the correct syntax would be?
Thanks!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.