Fits a shared parameter joint model for longitudinal and timetoevent (e.g. survival) data under a Bayesian framework using Stan.
stan_jm( formulaLong, dataLong, formulaEvent, dataEvent, time_var, id_var, family = gaussian, assoc = "etavalue", lag_assoc = 0, grp_assoc, epsilon = 1e05, basehaz = c("bs", "weibull", "piecewise"), basehaz_ops, qnodes = 15, init = "prefit", weights, priorLong = normal(autoscale = TRUE), priorLong_intercept = normal(autoscale = TRUE), priorLong_aux = cauchy(0, 5, autoscale = TRUE), priorEvent = normal(autoscale = TRUE), priorEvent_intercept = normal(autoscale = TRUE), priorEvent_aux = cauchy(autoscale = TRUE), priorEvent_assoc = normal(autoscale = TRUE), prior_covariance = lkj(autoscale = TRUE), prior_PD = FALSE, algorithm = c("sampling", "meanfield", "fullrank"), adapt_delta = NULL, max_treedepth = 10L, QR = FALSE, sparse = FALSE, ... )
formulaLong  A twosided linear formula object describing both the
fixedeffects and randomeffects parts of the longitudinal submodel,
similar in vein to formula specification in the lme4 package
(see 


dataLong  A data frame containing the variables specified in


formulaEvent  A twosided formula object describing the event
submodel. The left hand side of the formula should be a 

dataEvent  A data frame containing the variables specified in


time_var  A character string specifying the name of the variable
in 

id_var  A character string specifying the name of the variable in


family  The family (and possibly also the link function) for the
longitudinal submodel(s). See 

assoc  A character string or character vector specifying the joint
model association structure. Possible association structures that can
be used include: "etavalue" (the default); "etaslope"; "etaauc";
"muvalue"; "muslope"; "muauc"; "shared_b"; "shared_coef"; or "null".
These are described in the Details section below. For a multivariate
joint model, different association structures can optionally be used for
each longitudinal submodel by specifying a list of character
vectors, with each element of the list specifying the desired association
structure for one of the longitudinal submodels. Specifying 

lag_assoc  A nonnegative scalar specifying the time lag that should be
used for the association structure. That is, the hazard of the event at
time t will be assumed to be associated with the value/slope/auc of
the longitudinal marker at time tu, where u is the time lag.
If fitting a multivariate joint model, then a different time lag can be used
for each longitudinal marker by providing a numeric vector of lags, otherwise
if a scalar is provided then the specified time lag will be used for all
longitudinal markers. Note however that only one time lag can be specified
for linking each longitudinal marker to the
event, and that that time lag will be used for all association structure
types (e.g. 

grp_assoc  Character string specifying the method for combining information
across lower level units clustered within an individual when forming the
association structure. This is only relevant when a grouping factor is
specified in 

epsilon  The halfwidth of the central difference used to numerically
calculate the derivate when the 

basehaz  A character string indicating which baseline hazard to use
for the event submodel. Options are a Bsplines approximation estimated
for the log baseline hazard ( 

basehaz_ops  A named list specifying options related to the baseline
hazard. Currently this can include:


qnodes  The number of nodes to use for the GaussKronrod quadrature that is used to evaluate the cumulative hazard in the likelihood function. Options are 15 (the default), 11 or 7. 

init  The method for generating the initial values for the MCMC.
The default is 

weights  Experimental and should be used with caution. The user can optionally supply a 2column data frame containing a set of 'prior weights' to be used in the estimation process. The data frame should contain two columns: the first containing the IDs for each individual, and the second containing the corresponding weights. The data frame should only have one row for each individual; that is, weights should be constant within individuals. 

priorLong, priorEvent, priorEvent_assoc  The prior distributions for the regression coefficients in the longitudinal submodel(s), event submodel, and the association parameter(s). Can be a call to one of the various functions provided by rstanarm for specifying priors. The subset of these functions that can be used for the prior on the coefficients can be grouped into several "families":
See the priors help page for details on the families and
how to specify the arguments for all of the functions in the table above.
To omit a prior i.e., to use a flat (improper) uniform prior
Note: Unless 

priorLong_intercept, priorEvent_intercept  The prior distributions
for the intercepts in the longitudinal submodel(s) and event submodel.
Can be a call to Note: The prior distribution for the intercept is set so it applies to the value when all predictors are centered. Moreover, note that a prior is only placed on the intercept for the event submodel when a Weibull baseline hazard has been specified. For the Bsplines and piecewise constant baseline hazards there is not intercept parameter that is given a prior distribution; an intercept parameter will be shown in the output for the fitted model, but this just corresponds to the necessary postestimation adjustment in the linear predictor due to the centering of the predictiors in the event submodel. 

priorLong_aux  The prior distribution for the "auxiliary" parameters
in the longitudinal submodels (if applicable).
The "auxiliary" parameter refers to a different parameter
depending on the
If fitting a multivariate joint model, you have the option to specify a list of prior distributions, however the elements of the list that correspond to any longitudinal submodel which does not have an auxiliary parameter will be ignored. 

priorEvent_aux  The prior distribution for the "auxiliary" parameters
in the event submodel. The "auxiliary" parameters refers to different
parameters depending on the baseline hazard. For 

prior_covariance  Cannot be 

prior_PD  A logical scalar (defaulting to 

algorithm  A string (possibly abbreviated) indicating the
estimation approach to use. Can be 

adapt_delta  Only relevant if 

max_treedepth  A positive integer specifying the maximum treedepth
for the nonUturn sampler. See the 

QR  A logical scalar defaulting to 

sparse  A logical scalar (defaulting to 

...  Further arguments passed to the function in the rstan
package ( 
A stanjm object is returned.
The stan_jm
function can be used to fit a joint model (also
known as a shared parameter model) for longitudinal and timetoevent data
under a Bayesian framework. The underlying
estimation is carried out using the Bayesian C++ package Stan
(http://mcstan.org/).
The joint model may be univariate (with only one longitudinal submodel) or
multivariate (with more than one longitudinal submodel).
For the longitudinal submodel a (possibly multivariate) generalised linear
mixed model is assumed with any of the family
choices
allowed by glmer
. If a multivariate joint model is specified
(by providing a list of formulas in the formulaLong
argument), then
the multivariate longitudinal submodel consists of a multivariate generalized
linear model (GLM) with groupspecific terms that are assumed to be correlated
across the different GLM submodels. That is, within
a grouping factor (for example, patient ID) the groupspecific terms are
assumed to be correlated across the different GLM submodels. It is
possible to specify a different outcome type (for example a different
family and/or link function) for each of the GLM submodels, by providing
a list of family
objects in the family
argument. Multilevel
clustered data are allowed, and that additional clustering can occur at a
level higher than the individuallevel (e.g. patients clustered within
clinics), or at a level lower than the individuallevel (e.g. tumor lesions
clustered within patients). If the clustering occurs at a level lower than
the individual, then the user needs to indicate how the lower level
clusters should be handled when forming the association structure between
the longitudinal and event submodels (see the grp_assoc
argument
described above).
For the event submodel a parametric
proportional hazards model is assumed. The baseline hazard can be estimated
using either a cubic Bsplines approximation (basehaz = "bs"
, the
default), a Weibull distribution (basehaz = "weibull"
), or a
piecewise constant baseline hazard (basehaz = "piecewise"
).
If the Bspline or piecewise constant baseline hazards are used,
then the degrees of freedom or the internal knot locations can be
(optionally) specified. If
the degrees of freedom are specified (through the df
argument) then
the knot locations are automatically generated based on the
distribution of the observed event times (not including censoring times).
Otherwise internal knot locations can be specified
directly through the knots
argument. If neither df
or
knots
is specified, then the default is to set df
equal to 6.
It is not possible to specify both df
and knots
.
Timevarying covariates are allowed in both the
longitudinal and event submodels. These should be specified in the data
in the same way as they normally would when fitting a separate
longitudinal model using lmer
or a separate
timetoevent model using coxph
. These timevarying
covariates should be exogenous in nature, otherwise they would perhaps
be better specified as an additional outcome (i.e. by including them as an
additional longitudinal outcome in the joint model).
Bayesian estimation of the joint model is performed via MCMC. The Bayesian
model includes independent priors on the
regression coefficients for both the longitudinal and event submodels,
including the association parameter(s) (in much the same way as the
regression parameters in stan_glm
) and
priors on the terms of a decomposition of the covariance matrices of the
groupspecific parameters.
See priors
for more information about the priors distributions
that are available.
GaussKronrod quadrature is used to numerically evaluate the integral
over the cumulative hazard in the likelihood function for the event submodel.
The accuracy of the numerical approximation can be controlled using the
number of quadrature nodes, specified through the qnodes
argument. Using a higher number of quadrature nodes will result in a more
accurate approximation.
The association structure for the joint model can be based on any of the following parameterisations:
current value of the linear predictor in the
longitudinal submodel ("etavalue"
)
first derivative (slope) of the linear predictor in the
longitudinal submodel ("etaslope"
)
the area under the curve of the linear predictor in the
longitudinal submodel ("etaauc"
)
current expected value of the longitudinal submodel
("muvalue"
)
the area under the curve of the expected value from the
longitudinal submodel ("muauc"
)
shared individuallevel random effects ("shared_b"
)
shared individuallevel random effects which also incorporate
the corresponding fixed effect as well as any corresponding
random effects for clustering levels higher than the individual)
("shared_coef"
)
interactions between association terms and observed data/covariates
("etavalue_data"
, "etaslope_data"
, "muvalue_data"
,
"muslope_data"
). These are described further below.
interactions between association terms corresponding to different
longitudinal outcomes in a multivariate joint model
("etavalue_etavalue(#)"
, "etavalue_muvalue(#)"
,
"muvalue_etavalue(#)"
, "muvalue_muvalue(#)"
). These
are described further below.
no association structure (equivalent to fitting separate
longitudinal and event models) ("null"
or NULL
)
More than one association structure can be specified, however,
not all possible combinations are allowed.
Note that for the lagged association structures baseline values (time = 0)
are used for the instances
where the time lag results in a time prior to baseline. When using the
"etaauc"
or "muauc"
association structures, the area under
the curve is evaluated using GaussKronrod quadrature with 15 quadrature
nodes. By default, "shared_b"
and "shared_coef"
contribute
all random effects to the association structure; however, a subset of the
random effects can be chosen by specifying their indices between parentheses
as a suffix, for example, "shared_b(1)"
or "shared_b(1:3)"
or
"shared_b(1,2,4)"
, and so on.
In addition, several association terms ("etavalue"
, "etaslope"
,
"muvalue"
, "muslope"
) can be interacted with observed
data/covariates. To do this, use the association term's main handle plus a
suffix of "_data"
then followed by the model matrix formula in
parentheses. For example if we had a variable in our dataset for gender
named sex
then we might want to obtain different estimates for the
association between the current slope of the marker and the risk of the
event for each gender. To do this we would specify
assoc = c("etaslope", "etaslope_data(~ sex)")
.
It is also possible, when fitting a multivariate joint model, to include
interaction terms between the association terms themselves (this only
applies for interacting "etavalue"
or "muvalue"
). For example,
if we had a joint model with two longitudinal markers, we could specify
assoc = list(c("etavalue", "etavalue_etavalue(2)"), "etavalue")
.
The first element of list says we want to use the value of the linear
predictor for the first marker, as well as it's interaction with the
value of the linear predictor for the second marker. The second element of
the list says we want to also include the expected value of the second marker
(i.e. as a "main effect"). Therefore, the linear predictor for the event
submodel would include the "main effects" for each marker as well as their
interaction.
There are additional examples in the Examples section below.
stanregobjects
, stanmvregmethods
,
print.stanmvreg
, summary.stanmvreg
,
posterior_traj
, posterior_survfit
,
posterior_predict
, posterior_interval
,
pp_check
, ps_check
, stan_mvmer
.
# \donttest{ if (.Platform$OS.type != "windows"  .Platform$r_arch !="i386") { ##### # Univariate joint model, with association structure based on the # current value of the linear predictor f1 < stan_jm(formulaLong = logBili ~ year + (1  id), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, time_var = "year", # this next line is only to keep the example small in size! chains = 1, cores = 1, seed = 12345, iter = 1000) print(f1) summary(f1) ##### # Univariate joint model, with association structure based on the # current value and slope of the linear predictor f2 < stan_jm(formulaLong = logBili ~ year + (year  id), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, assoc = c("etavalue", "etaslope"), time_var = "year", chains = 1, cores = 1, seed = 12345, iter = 1000) print(f2) ##### # Univariate joint model, with association structure based on the # lagged value of the linear predictor, where the lag is 2 time # units (i.e. 2 years in this example) f3 < stan_jm(formulaLong = logBili ~ year + (1  id), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, time_var = "year", assoc = "etavalue", lag_assoc = 2, chains = 1, cores = 1, seed = 12345, iter = 1000) print(f3) ##### # Univariate joint model, where the association structure includes # interactions with observed data. Here we specify that we want to use # an association structure based on the current value of the linear # predictor from the longitudinal submodel (i.e. "etavalue"), but we # also want to interact this with the treatment covariate (trt) from # pbcLong data frame, so that we can estimate a different association # parameter (i.e. estimated effect of log serum bilirubin on the log # hazard of death) for each treatment group f4 < stan_jm(formulaLong = logBili ~ year + (1  id), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, time_var = "year", assoc = c("etavalue", "etavalue_data(~ trt)"), chains = 1, cores = 1, seed = 12345, iter = 1000) print(f4) ###### # Multivariate joint model, with association structure based # on the current value and slope of the linear predictor in the # first longitudinal submodel and the area under the marker # trajectory for the second longitudinal submodel mv1 < stan_jm( formulaLong = list( logBili ~ year + (1  id), albumin ~ sex + year + (year  id)), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, assoc = list(c("etavalue", "etaslope"), "etaauc"), time_var = "year", chains = 1, cores = 1, seed = 12345, iter = 100) print(mv1) ##### # Multivariate joint model, where the association structure is formed by # including the expected value of each longitudinal marker (logBili and # albumin) in the linear predictor of the event submodel, as well as their # interaction effect (i.e. the interaction between the two "etavalue" terms). # Note that whether such an association structure based on a marker by # marker interaction term makes sense will depend on the context of your # application  here we just show it for demostration purposes). mv2 < stan_jm( formulaLong = list( logBili ~ year + (1  id), albumin ~ sex + year + (year  id)), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, assoc = list(c("etavalue", "etavalue_etavalue(2)"), "etavalue"), time_var = "year", chains = 1, cores = 1, seed = 12345, iter = 100) ##### # Multivariate joint model, with one bernoulli marker and one # Gaussian marker. We will artificially create the bernoulli # marker by dichotomising log serum bilirubin pbcLong$ybern < as.integer(pbcLong$logBili >= mean(pbcLong$logBili)) mv3 < stan_jm( formulaLong = list( ybern ~ year + (1  id), albumin ~ sex + year + (year  id)), dataLong = pbcLong, formulaEvent = Surv(futimeYears, death) ~ sex + trt, dataEvent = pbcSurv, family = list(binomial, gaussian), time_var = "year", chains = 1, cores = 1, seed = 12345, iter = 1000) }#> Fitting a univariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000255 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.55 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 4.85136 seconds (Warmup) #> Chain 1: 3.04527 seconds (Sampling) #> Chain 1: 7.89663 seconds (Total) #> Chain 1:#> Warning: The largest Rhat is 1.05, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#rhat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailess#> stan_jm #> formula (Long1): logBili ~ year + (1  id) #> family (Long1): gaussian [identity] #> formula (Event): Surv(futimeYears, death) ~ sex + trt #> baseline hazard: bs #> assoc: etavalue (Long1) #>  #> #> Longitudinal submodel: logBili #> Median MAD_SD #> (Intercept) 0.838 0.213 #> year 0.092 0.010 #> sigma 0.509 0.021 #> #> Event submodel: #> Median MAD_SD exp(Median) #> (Intercept) 2.991 0.649 0.050 #> sexf 0.397 0.609 0.672 #> trt 0.746 0.477 0.474 #> Long1etavalue 1.425 0.269 4.158 #> bsplinescoef1 1.238 0.975 NA #> bsplinescoef2 0.087 0.885 NA #> bsplinescoef3 1.544 1.270 NA #> bsplinescoef4 0.513 1.510 NA #> bsplinescoef5 0.097 1.768 NA #> bsplinescoef6 0.457 1.657 NA #> #> Grouplevel error terms: #> Groups Name Std.Dev. #> id Long1(Intercept) 1.364 #> Num. levels: id 40 #> #> Sample avg. posterior predictive distribution #> of longitudinal outcomes: #> Median MAD_SD #> Long1mean_PPD 0.584 0.042 #> #>  #> For info on the priors used see help('prior_summary.stanreg').Fitting a univariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000346 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.46 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 26.3729 seconds (Warmup) #> Chain 1: 16.4721 seconds (Sampling) #> Chain 1: 42.845 seconds (Total) #> Chain 1:#> Warning: The largest Rhat is 1.05, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#rhat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> stan_jm #> formula (Long1): logBili ~ year + (year  id) #> family (Long1): gaussian [identity] #> formula (Event): Surv(futimeYears, death) ~ sex + trt #> baseline hazard: bs #> assoc: etavalue (Long1), etaslope (Long1) #>  #> #> Longitudinal submodel: logBili #> Median MAD_SD #> (Intercept) 0.683 0.190 #> year 0.272 0.053 #> sigma 0.356 0.016 #> #> Event submodel: #> Median MAD_SD exp(Median) #> (Intercept) 3.407 0.818 0.033 #> sexf 0.383 0.675 0.682 #> trt 1.015 0.563 0.362 #> Long1etavalue 0.771 0.427 2.163 #> Long1etaslope 10.463 5.376 35012.688 #> bsplinescoef1 6.099 4.257 NA #> bsplinescoef2 2.506 2.159 NA #> bsplinescoef3 3.166 1.613 NA #> bsplinescoef4 0.316 1.740 NA #> bsplinescoef5 0.218 1.757 NA #> bsplinescoef6 0.527 1.696 NA #> #> Grouplevel error terms: #> Groups Name Std.Dev. Corr #> id Long1(Intercept) 1.2680 #> Long1year 0.2507 0.69 #> Num. levels: id 40 #> #> Sample avg. posterior predictive distribution #> of longitudinal outcomes: #> Median MAD_SD #> Long1mean_PPD 0.586 0.028 #> #>  #> For info on the priors used see help('prior_summary.stanreg').Fitting a univariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000192 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.92 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 4.21702 seconds (Warmup) #> Chain 1: 2.5458 seconds (Sampling) #> Chain 1: 6.76282 seconds (Total) #> Chain 1:#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailess#> stan_jm #> formula (Long1): logBili ~ year + (1  id) #> family (Long1): gaussian [identity] #> formula (Event): Surv(futimeYears, death) ~ sex + trt #> baseline hazard: bs #> assoc: etavalue (Long1) #>  #> #> Longitudinal submodel: logBili #> Median MAD_SD #> (Intercept) 0.772 0.246 #> year 0.091 0.010 #> sigma 0.508 0.022 #> #> Event submodel: #> Median MAD_SD exp(Median) #> (Intercept) 2.974 0.585 0.051 #> sexf 0.312 0.589 0.732 #> trt 0.651 0.423 0.521 #> Long1etavalue 1.432 0.267 4.185 #> bsplinescoef1 1.474 1.132 NA #> bsplinescoef2 0.174 0.881 NA #> bsplinescoef3 1.298 1.121 NA #> bsplinescoef4 0.591 1.493 NA #> bsplinescoef5 0.183 1.664 NA #> bsplinescoef6 0.432 1.818 NA #> #> Grouplevel error terms: #> Groups Name Std.Dev. #> id Long1(Intercept) 1.319 #> Num. levels: id 40 #> #> Sample avg. posterior predictive distribution #> of longitudinal outcomes: #> Median MAD_SD #> Long1mean_PPD 0.585 0.041 #> #>  #> For info on the priors used see help('prior_summary.stanreg').Fitting a univariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000227 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.27 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 6.07398 seconds (Warmup) #> Chain 1: 3.02875 seconds (Sampling) #> Chain 1: 9.10273 seconds (Total) #> Chain 1:#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailess#> stan_jm #> formula (Long1): logBili ~ year + (1  id) #> family (Long1): gaussian [identity] #> formula (Event): Surv(futimeYears, death) ~ sex + trt #> baseline hazard: bs #> assoc: etavalue (Long1), etavalue_data (Long1) #>  #> #> Longitudinal submodel: logBili #> Median MAD_SD #> (Intercept) 0.793 0.222 #> year 0.092 0.010 #> sigma 0.507 0.025 #> #> Event submodel: #> Median MAD_SD exp(Median) #> (Intercept) 3.144 0.607 0.043 #> sexf 0.391 0.555 0.676 #> trt 0.329 0.908 0.720 #> Long1etavalue 1.540 0.345 4.667 #> Long1etavalue:trt 0.224 0.562 0.799 #> bsplinescoef1 1.407 1.081 NA #> bsplinescoef2 0.111 0.803 NA #> bsplinescoef3 1.725 1.129 NA #> bsplinescoef4 0.754 1.491 NA #> bsplinescoef5 0.221 1.498 NA #> bsplinescoef6 0.282 1.595 NA #> #> Grouplevel error terms: #> Groups Name Std.Dev. #> id Long1(Intercept) 1.298 #> Num. levels: id 40 #> #> Sample avg. posterior predictive distribution #> of longitudinal outcomes: #> Median MAD_SD #> Long1mean_PPD 0.586 0.040 #> #>  #> For info on the priors used see help('prior_summary.stanreg').Fitting a multivariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: Rejecting initial value: #> Chain 1: Log probability evaluates to log(0), i.e. negative infinity. #> Chain 1: Stan can't start sampling from this initial value. #> Chain 1: #> Chain 1: Gradient evaluation took 0.003503 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 35.03 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: WARNING: There aren't enough warmup iterations to fit the #> Chain 1: three stages of adaptation as currently configured. #> Chain 1: Reducing each adaptation stage to 15%/75%/10% of #> Chain 1: the given number of warmup iterations: #> Chain 1: init_buffer = 7 #> Chain 1: adapt_window = 38 #> Chain 1: term_buffer = 5 #> Chain 1: #> Chain 1: Iteration: 1 / 100 [ 1%] (Warmup) #> Chain 1: Iteration: 10 / 100 [ 10%] (Warmup) #> Chain 1: Iteration: 20 / 100 [ 20%] (Warmup) #> Chain 1: Iteration: 30 / 100 [ 30%] (Warmup) #> Chain 1: Iteration: 40 / 100 [ 40%] (Warmup) #> Chain 1: Iteration: 50 / 100 [ 50%] (Warmup) #> Chain 1: Iteration: 51 / 100 [ 51%] (Sampling) #> Chain 1: Iteration: 60 / 100 [ 60%] (Sampling) #> Chain 1: Iteration: 70 / 100 [ 70%] (Sampling) #> Chain 1: Iteration: 80 / 100 [ 80%] (Sampling) #> Chain 1: Iteration: 90 / 100 [ 90%] (Sampling) #> Chain 1: Iteration: 100 / 100 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 39.8348 seconds (Warmup) #> Chain 1: 43.8592 seconds (Sampling) #> Chain 1: 83.694 seconds (Total) #> Chain 1:#> Warning: There were 30 divergent transitions after warmup. Increasing adapt_delta above 0.85 may help. See #> http://mcstan.org/misc/warnings.html#divergenttransitionsafterwarmup#> Warning: There were 20 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See #> http://mcstan.org/misc/warnings.html#maximumtreedepthexceeded#> Warning: Examine the pairs() plot to diagnose sampling problems#> Warning: The largest Rhat is 2.12, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#rhat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailess#> Warning: Markov chains did not converge! Do not analyze results!#> stan_jm #> formula (Long1): logBili ~ year + (1  id) #> family (Long1): gaussian [identity] #> formula (Long2): albumin ~ sex + year + (year  id) #> family (Long2): gaussian [identity] #> formula (Event): Surv(futimeYears, death) ~ sex + trt #> baseline hazard: bs #> assoc: etavalue (Long1), etaslope (Long1), etaauc (Long2) #>  #> #> Longitudinal submodel 1: logBili #> Median MAD_SD #> (Intercept) 0.697 0.000 #> year 0.086 0.000 #> sigma 0.526 0.000 #> #> Longitudinal submodel 2: albumin #> Median MAD_SD #> (Intercept) 3.509 0.000 #> sexf 0.081 0.000 #> year 0.123 0.000 #> sigma 0.334 0.000 #> #> Event submodel: #> Median MAD_SD exp(Median) #> (Intercept) 1.082885e+11 5.004000e+01 Inf #> sexf 3.470000e01 0.000000e+00 7.070000e01 #> trt 9.600000e02 0.000000e+00 9.080000e01 #> Long1etavalue 3.422000e+00 0.000000e+00 3.064300e+01 #> Long1etaslope 1.264015e+12 5.840970e+02 0.000000e+00 #> Long2etaauc 3.360000e01 0.000000e+00 1.399000e+00 #> bsplinescoef1 0.000000e+00 0.000000e+00 NA #> bsplinescoef2 0.000000e+00 0.000000e+00 NA #> bsplinescoef3 0.000000e+00 0.000000e+00 NA #> bsplinescoef4 0.000000e+00 0.000000e+00 NA #> bsplinescoef5 0.000000e+00 0.000000e+00 NA #> bsplinescoef6 0.000000e+00 0.000000e+00 NA #> #> Grouplevel error terms: #> Groups Name Std.Dev. Corr #> id Long1(Intercept) 1.11688 #> Long2(Intercept) 0.36463 0.00 #> Long2year 0.03419 0.00 0.00 #> Num. levels: id 40 #> #> Sample avg. posterior predictive distribution #> of longitudinal outcomes: #> Median MAD_SD #> Long1mean_PPD 0.878 0.035 #> Long2mean_PPD 3.142 0.014 #> #>  #> For info on the priors used see help('prior_summary.stanreg').Fitting a multivariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000441 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.41 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: WARNING: There aren't enough warmup iterations to fit the #> Chain 1: three stages of adaptation as currently configured. #> Chain 1: Reducing each adaptation stage to 15%/75%/10% of #> Chain 1: the given number of warmup iterations: #> Chain 1: init_buffer = 7 #> Chain 1: adapt_window = 38 #> Chain 1: term_buffer = 5 #> Chain 1: #> Chain 1: Iteration: 1 / 100 [ 1%] (Warmup) #> Chain 1: Iteration: 10 / 100 [ 10%] (Warmup) #> Chain 1: Iteration: 20 / 100 [ 20%] (Warmup) #> Chain 1: Iteration: 30 / 100 [ 30%] (Warmup) #> Chain 1: Iteration: 40 / 100 [ 40%] (Warmup) #> Chain 1: Iteration: 50 / 100 [ 50%] (Warmup) #> Chain 1: Iteration: 51 / 100 [ 51%] (Sampling) #> Chain 1: Iteration: 60 / 100 [ 60%] (Sampling) #> Chain 1: Iteration: 70 / 100 [ 70%] (Sampling) #> Chain 1: Iteration: 80 / 100 [ 80%] (Sampling) #> Chain 1: Iteration: 90 / 100 [ 90%] (Sampling) #> Chain 1: Iteration: 100 / 100 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 5.3195 seconds (Warmup) #> Chain 1: 14.0857 seconds (Sampling) #> Chain 1: 19.4052 seconds (Total) #> Chain 1:#> Warning: There were 7 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See #> http://mcstan.org/misc/warnings.html#maximumtreedepthexceeded#> Warning: Examine the pairs() plot to diagnose sampling problems#> Warning: The largest Rhat is 1.3, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#rhat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#tailess#> Warning: Markov chains did not converge! Do not analyze results!#> Fitting a multivariate joint model. #> #> Please note the warmup may be much slower than later iterations! #> #> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.000351 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.51 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 23.3267 seconds (Warmup) #> Chain 1: 10.9445 seconds (Sampling) #> Chain 1: 34.2712 seconds (Total) #> Chain 1:#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mcstan.org/misc/warnings.html#bulkess# }