Default baseline hazard knots format and placement when Bsplines_degree > 0L

[edbonneville]

Hi {JMbayes2} team - thank you for the impressive work on this package!

Related to #24, but I was noticing strange results (and totally different compared to {JM}) when using the baseline hazards that were not piecewise exponential. It was with an applied dataset where events only started occurring after 1 month (and maximum follow-up was 6 months, with administrative censoring). I tried to reproduce a similar set up based on the pbc2 data below:

set.seed(20220715)
library(JM)
#> Loading required package: MASS
#> Loading required package: nlme
#> Loading required package: splines
#> Loading required package: survival
library(JMbayes2)
#> Loading required package: GLMMadaptive
#> 
#> Attaching package: 'GLMMadaptive'
#> The following object is masked from 'package:MASS':
#> 
#>     negative.binomial
#> 
#> Attaching package: 'JMbayes2'
#> The following object is masked from 'package:MASS':
#> 
#>     area

# Suppose events only start after 3 years
ids_3y <- with(pbc2.id, id[which(years >= 3)])
pbc2.id_3y <- subset(pbc2.id, id %in% ids_3y)
pbc2_3y <- subset(pbc2, id %in% ids_3y)

# Fit submodels - second example from JM::jointModel()
lmeFit <- lme(log(serBilir) ~ year * drug, random = ~ year | id, data = pbc2_3y)
coxFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id_3y, x = TRUE)

# Fit with JM
pbc_JM <- jointModel(
  lmeObject = lmeFit,
  survObject = coxFit,
  method = "spline-PH-aGH",
  timeVar = "year",
  ord = 3L, # quadratic B-splines
  lng.in.kn = 9L # internal knots, so 10 "base_hazard_segments"
)

# Check knots (boundary knots repeated ord = 3L times)
pbc_JM$control$knots
#> $`1`
#>  [1]  0.01347527  0.01347527  0.01347527  4.21202497  5.18946446  5.73978754
#>  [7]  6.51626328  7.21443434  8.19789727  8.98368196 10.03121235 12.16679444
#> [13] 14.30566203 14.30566203 14.30566203

# Fit JMbayes2 with defaults
pbc_JMbayes2_m1 <- jm(
  Mixed_objects = lmeFit,
  Surv_object = coxFit,
  time_var = "year",
  Bsplines_degree = 2L, # quadratic, the default
  base_hazard_segments = 10L # also default
)

# Check knots
pbc_JMbayes2_m1$control$knots
#> $`1`
#>  [1]  0.9237761  2.0389333  3.1540905  4.2692476  5.3844048  6.4995619
#>  [7]  7.6147191  8.7298762  9.8450334 10.9601906 12.0753477 13.1905049
#> [13] 14.3056620 15.4208192 16.5359763

# Fit JMbayes2 with knots as formatted by JM
pbc_JMbayes2_m2 <- jm(
  Mixed_objects = lmeFit,
  Surv_object = coxFit,
  time_var = "year",
  knots = pbc_JM$control$knots,
  Bsplines_degree = 2L
)

# Compare coefficients
cbind.data.frame(
  "JM" = coef(pbc_JM, "Event"),
  "JMbayes2 (own knots)" = c(
    coef(pbc_JMbayes2_m1)$gammas,
    coef(pbc_JMbayes2_m1)$association
  ),
  "JMbayes2 (JM knots)" = c(
    coef(pbc_JMbayes2_m2)$gammas,
    coef(pbc_JMbayes2_m2)$association
  )
)
#>                       JM JMbayes2 (own knots) JMbayes2 (JM knots)
#> drugD-penicil 0.09590339            0.0196705          0.07466035
#> Assoct        1.09028319            0.6905448          1.12973032

^{Created on 2022-07-15 by the reprex package (v2.0.1)}

So beyond the results being very different, there is also:

1. Some of the knots from pbc_JMbayes2_m1 are placed outside of the follow-up window. If you try the same code with Bsplines_degree = 3L, the first knot will even be placed at a negative timepoint!
2. The format of the knots (I think created here) does not seem to be the one required when passed to splineDesign() here. Based on what JM::jointModel() does, it would seem you need the boundary knots to be repeated Bsplines_degree + 1L times. Somehow though, things still work out with the defaults if you have data with events starting from the beginning of follow-up (if you re-run the above example with the original pbc2 datasets, results are the same for all three models..)

I am using the latest version of {JM}, and the latest development version of {JMbayes2}. Thanks in advance!

[drizopoulos]

Thanks for reporting this. The issues should be resolved in the development version on GitHub. We will do some extra tests before releasing it to CRAN.

[edbonneville]

Thanks a lot! I have not tried it yet but I wonder how the development version will now choose between two identically named knots() helper functions.. no problem waiting for the next CRAN release.