drizopoulos / jm

Joint Models for Longitudinal & Survival Data under Maximum Likelihood

JM: Joint Models for Longitudinal and Survival Data using Maximum Likelihood

Description

This repository contains the source files for the R package JM. This package fits joint models for longitudinal and time-to-event data using maximum likelihood. These models are applicable in mainly two settings. First, when focus is on the survival outcome and we wish to account for the effect of an endogenous (aka internal) time-dependent covariates measured with error. Second, when focus is on the longitudinal outcome and we wish to correct for nonrandom dropout.

The basic joint-model-fitting function of the package is jointModel(). This accepts as main arguments a linear mixed model fitted by function lme() from the nlme package and a Cox model fitted using function coxph() from the survival package.

Basic Features

• It can fit joint models for a single continuous longitudinal outcome and a time-to-event outcome.

• For the survival outcome a relative risk models is assumed. The method argument of jointModel() can be used to define the type of baseline hazard function. Options are a B-spline approximation, a piecewise-constant function, the Weibull hazard and a completely unspecified function (i.e., a discrete function with point masses at the unique event times).

• The user has now the option to define custom transformation functions for the terms of the longitudinal submodel that enter into the linear predictor of the survival submodel (arguments derivForm, parameterization). For example, the current value of the longitudinal outcomes, the velocity of the longitudinal outcome (slope), the area under the longitudinal profile. From the aforementioned options, in each model up to two terms can be included. In addition, using argument InterFact interactions terms can be considered.

Dynamic predictions

• Function survfitJM() computes dynamic survival probabilities.

• Function predict() computes dynamic predictions for the longitudinal outcome.

• Function aucJM() calculates time-dependent AUCs for joint models, and function rocJM() calculates the corresponding time-dependent sensitivities and specifies.

• Function prederrJM() calculates prediction errors for joint models.

jm's Issues

Confusing point estimates for association parameters in current values + current slopes model

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!

hard coded n of knots is causing an error with custom knots

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:

JM/R/log.posterior.b.R

Line 62 in 05bfde1

wkP <- wk * rep(P, each = 7)

problem running fuction jointModel() in JM package

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)
summary(model.cox)
fit.lm<-lme(logTSH time,random= ~time | code,data=dataLong)
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,

`survfitJM` expects newdata to be ordered by subject, and by time within each subject, else things go wrong

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.

Warning on second run of predict.jointModel

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.)

jointModel() looks for external data.frame

... 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...

problem running fuction jointModel() in JM package with method = "Cox-PH-aGH" option

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)

the linear mixed model

lmeFit1 <- lme(lserBilir ~ year + drug:year, data = pbc2,
random = ~ year | id)

the Cox model

survFit1 <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

the joint model

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!

Residuals for competing risks models

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.

Infinite or missing values in Hessian at convergence.

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")

Error in exp(coefs[1]) : non-numeric argument to mathematical function- fail to run CompRisk

I try to run sample codes from http://jmr.r-forge.r-project.org/Chapter5.html

Section 5.5.1

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!

Error in if (w_dif < specs$tol || iter == specs$maxiter) break : missing value where TRUE/FALSE needed

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

Trouble with specifying integral form for cumulative effects model

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),
random=-1+time+time^2/2,indRandom=c(1:2))>

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!

problem with scaling of time variables

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?

rocJM producing unexpected effects

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.

Error in aeqSurv(Y) : aeqSurv exception, an interval has effective length 0

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")

error when applying weight function to ins()

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))+
ins(time,knots=3,weight.fun =g)+
I(g(time)*(cov)),
indFixed=1:5,
random=-1+I(g(time))+I(g(time,1)),indRandom=1:2)

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!

***caught segfault*** caused by Error in names(fit$coef) <- cname : 'names' attribute [4] must be the same length as the vector [2]

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

Potential bug in survfitJM.jointModel.R

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

Error in if (t1 || t2) { : missing value where TRUE/FALSE needed

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!

“aeqSurv exception, an interval has effective length 0”

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!

survfitJM error: Error in rowMeans(rr, na.rm = TRUE) : 'x' must be numeric

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

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:

Simulating data and GLMMPQL using JM

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

Error in aeqSurv(Y)

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])

the cluster slot in survobj has been renamed to '(cluster)'

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)

Error in in optim(thetas, LogLik.splineGH, Score.splineGH, method = "BFGS",

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.

Error in solve.default(VC) : system is computationally singular: reciprocal condition number = 3.40196e-17

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.

### Please, if someone could help me;; i'm really stuck

Linear mixed model

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)

Cox

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

JOINT MODEL

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)

diagnostic plots

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!

Can the joint modeling frame work accept patients with different start time ?

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.

issues in running survfitJM

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

Fit a joint model with competing risks data

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

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

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

plot.survfitJM in shiny app : comparison of these types is not implemented

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

Extract fitted lme object from the jointModel

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?

Error in FUN(X[[i]], ...) : subscript out of bounds

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)

diagnostic plots

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 :)

Error while iterating: missing value where TRUE/FALSE needed

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

code graphical time dependent slope

The code for Figure 5-1 in the book "Joint Models for Longitudinal and Survival Data" on page 103 is as follows

Error in optim(thetas, opt.survPC, gr.survPC, method = "BFGS", control = list(maxit = if (it < : non-finite value supplied by optim

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.

Error in aeqSurv(Y) : aeqSurv exception, an interval has effective length 0

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!

Query on using `survfitJM`

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.

sample sizes in the longitudinal and event processes differ.

          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)

Query on getting back association parameter from data simulation in weibull-ph-agh

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

Question about the calculation of standard deviation in the "random effects parameterization"

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!

The interpretation of exp(Assoct) in JM

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

Error in Xderiv %*% fixef(lmeObject)[derivForm$indFixed] : non-conformable arguments

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

Accept models from lme4 in addition to nlme

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!

Comparison of floating point numbers

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).