Using random slopes but with fixed intercept

Hi! I was wondering what happens behind the scenes when running a model with random slopes but fixed intercept? The results between {JM} and {JMbayes2} are quite different in the example below using the pbc2 data:

set.seed(202207151)
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

# Fit submodel with random slopes but fixed intercept
lmeFit <- lme(log(serBilir) ~ ns(year, 2) * drug, random = ~ 0 + ns(year, 2) | id, data = pbc2)
coxFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

# JM fit
pbc_JM <- jointModel(
  lmeObject = lmeFit,
  survObject = coxFit,
  method = "piecewise-PH-aGH",
  timeVar = "year",
  lng.in.kn = 10L
)

# JMbayes2 fit
pbc_JMbayes2 <- jm(
  Mixed_objects = lmeFit,
  Surv_object = coxFit,
  time_var = "year",
  Bsplines_degree = 0L,
  knots = list(c(0, pbc_JM$control$knots, max(pbc2.id$years) + 1))
)

# Compare coefficients
cbind.data.frame(
  "JM" = coef(pbc_JM, "Event"),
  "JMbayes2" = c(
    coef(pbc_JMbayes2)$gammas,
    coef(pbc_JMbayes2)$association
  )
)
#>                       JM   JMbayes2
#> drugD-penicil 0.01240778 -0.0389509
#> Assoct        1.34275739  0.7221391

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

The coefficients for the survival submodels are different, but so to are the variance-covariance matrices of the random effects (despite being on same dimensions). I also used splines the example above because using linear random slopes and fixed intercept yields an error:

# Using just linear random slope
update(
  object = pbc_JMbayes2,
  Mixed_objects = update(lmeFit, fixed. = . ~ year * drug, random = ~ 0 + year | id)
)
#> Error in jm(Surv_object = coxFit, Mixed_objects = update(lmeFit, fixed. = . ~ : jm() does not currently work when you have a single longitudinal outcome and only random intercepts.

Note I get the same results as model pbc_JM using rstanarm::stan_jm() as below:

library(rstanarm)
options(mc.cores = 3)

pbc_stan <- stan_jm(
  formulaLong = log(serBilir) ~ ns(year, 2) * drug + (0 + ns(year, 2) | id),
  dataLong = pbc2,
  formulaEvent = survival::Surv(years, status2) ~ drug,
  dataEvent = pbc2.id,
  time_var = "year",
  basehaz = "piecewise",
  basehaz_ops = list("knots" = pbc_JM$control$knots),
  chains = 3,
  cores = 3
)

Thanks again in advance!

drizopoulos / JMbayes2

Using random slopes but with fixed intercept #32