library("JMbayes2")

pbc2.idCR <- crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", 
                         nameStrata = "CR")
CoxFit_CR <- coxph(Surv(years, status2) ~ (age + drug) * strata(CR),
                   data = pbc2.idCR)

fm1 <- lme(log(serBilir) ~ poly(year, 2) * drug, data = pbc2, 
           random = ~ poly(year, 2) | id)
fm2 <- lme(prothrombin ~ year * drug, data = pbc2, random = ~ year | id)

CR_forms <- list(
  "log(serBilir)" = ~ value(log(serBilir)):CR,
  "prothrombin" = ~ value(prothrombin):CR
)

# A short-run model for the purpose of generating default priors
jFit_CR <- jm(CoxFit_CR, list(fm1, fm2), time_var = "year", 
              functional_forms = CR_forms, 
              n_iter = 250L, n_burnin = 200L, n_thin = 1L)

# Set Tau_prior on 'year' in fm2 ...
#   ... you get much slower convergence and wildly biased estimates...
#   ... but eventually, the standard error is suitable small
# The narrow prior seems to really throw off the MCMC sampling, and I
#   would argue that the parameter space is not being searched well at all
Tau_prior <- jFit_CR$priors$Tau_betas_HC
largeTau  <- Tau_prior %*% diag(c(1,1,1,1,1,1,1,1e3,1,1))

# Note: this model takes much longer to converge
jFit_CR_largeTau <- jm(CoxFit_CR, list(fm1, fm2), time_var = "year", 
              functional_forms = CR_forms, 
              n_iter = 25000L, n_burnin = 5000L, n_thin = 5L,
              priors = list("Tau_betas_HC" = largeTau)
)

summary(jFit_CR_largeTau)$Outcome1
summary(jFit_CR_largeTau)$Outcome2

# Note: initial values are same for both models
# For largeTau, posterior means are very large but do not explore the parameter space much at all

par(mfrow = c(5,2))
traceplot(jFit_CR_largeTau, parm = "betas")