McFetridge et al

*Özgür Asar, Department of Biostatistics and Medical Informatics, Faculty of Medicine, Acıbadem Mehmet Ali Aydınlar University, İstanbul, Turkey.

# Robust joint modelling of longitudinal and survival data with a time-varying
degrees-of-freedom parameter

Lisa M. McFetridge$\dagger$ Özgür Asar*$\dagger$ Jonas Wallin  Mathematical Sciences Research Centre, School of Mathematics and Physics, Queen’s University of Belfast, Belfast, United Kingdom  Department of Biostatistics and Medical Informatics, Faculty of Medicine, Acıbadem Mehmet Ali Aydınlar University, İstanbul, Turkey  Department of Statistics, Lund University, Lund, Sweden<EMAIL_ADDRESS>|<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary]Repeated measures of biomarkers have the potential of explaining
hazards of survival outcomes. In practice, these measurements are
intermittently measured and are known to be subject to substantial measurement
error. Joint modelling of longitudinal and survival data enables us to
associate intermittently measured error-prone biomarkers with risks of
survival outcomes. Most of the joint models available in the literature have
been built on the Gaussian assumption. This makes them sensitive to outliers.
In this work, we study a range of robust models to address this issue. For
medical data, it has been observed that outliers might occur with different
frequencies over time. To address this, a new model with a time varying
robustness is introduced. Through both a simulation study and analysis of two
real-life data examples, this research not only stresses the need to account
for longitudinal outliers in joint modelling research but also highlights the
bias and inefficiency from not properly estimating the degrees-of-freedom
parameter. Each technique presented in this work can be fitted using the R
package robjm.

###### keywords:

Degrees-of-freedom; longitudinal outliers; normal variance-mixtures; robust
joint model; t-distribution

††articletype: Article Type00footnotetext: † L.M. McFetridge and Ö. Asar are
equal contributing authors.

##  1 Introduction

Biomarkers are used as proxies of one’s health, and are known to be subject to
substantial measurement error, due to biological and non-biological sources.
In prospective studies, the biomarkers are repeatedly measured. These
measurements are made at intermittent time-points in practice rather than
continuously, and times elapsed between successive measurements are typically
unequal. It is of scientific interest to explain risks of survival outcomes by
repeated measures of biomarker data. Two main obstacles are the intermittent
nature of the biomarker data and inherent measurement error. Joint modelling
of longitudinal and survival outcomes enables us to associate repeated
measures of biomarkers with the risks of survival events, while taking into
account the aforementioned obstacles. The framework typically combines two
sub-models, a mixed-effects model for repeated measures, and a Cox model for
survival data. The sub-models are linked with shared parameters. Reviews of
relevant literature can be found in Tsiatis and Davidian 1, Diggle, Sousa and
Chetwynd 2, Rizopoulos 3, McCrink, Marshall and Cairns 4, Asar et al. 5 and
Elashoff 6.

A prevailing assumption in the literature is that both random-effects and
error terms in the mixed-effects sub-model follow Normal distributions.
However, in most real-life problems, the data-sets include outliers for which
the Normal assumption might be inadequate. Specifically, in prospective
studies, two types of outliers may be present:

  1. i)

b-outliers: Outlying individuals within the population that do not conform to
population trends. These are outliers within the longitudinal random-effects.

  2. ii)

e-outliers: Outlying observations within an individual’s set of measurements
that do not follow the individual’s own trend over time. These are outliers
within the longitudinal random error.

Only in recent years have studies been undertaken to investigate the negative
impacts that the presence of longitudinal outliers under the Normal assumption
can have 7, 8, 9, 10, 11. Robust joint models replace the Normality assumption
with t-distributional assumptions for the random terms. In doing so, the
heavier tails down-weigh the detrimental impact of longitudinal outliers
within the joint modelling framework. Initial research on robust joint
modelling considered only the $e$-outliers 7, 8, 12, and both $e$\- and
$b$-outliers 10, 11, 13, 14. Each of these studies found that the parameter
estimates and corresponding standard errors are sensitive to the Normality
assumptions when outliers are present, with t-distributional assumptions
alleviating such bias. If utilised when outliers are not present, unbiased
estimates are obtained, though Li, Elashoff and Li 7 and Huang, Li and
Elashoff 8 noted that the robust joint models give slightly higher standard
errors for the longitudinal parameters. This initial work however
restrictively fixes the degrees-of-freedom at a constant chosen by the user,
an assumption alleviated by Baghfalaki, Ganjali and Berridge 9, Asar, Fournier
and Dantan 11 and Baghfalaki, Ganjali and Hashemi 14, who utilised Bayesian
approaches, and McCrink 10, who utilised a frequentist approach, to allow
estimation of the degrees-of-freedom to be dictated by the data. With
estimated degrees-of-freedom, the t-distributional assumptions of the robust
joint model will approximate normality in the absence of outliers.

All of the aforementioned research, however, assumes that the degrees-of-
freedom is constant, unchanging over time, an assumption prevalent in the
literature for robust mixed models alongside their robust joint model
counterparts. To address this issue, we propose a robust joint model which can
account for the situation where the impact of outliers can change over time.
Varying frequency of outliers is a very likely scenario, for example, when
patients are given new treatments, they typically take time to stabilise and
adjust to these treatments. This period of adjustment means that they are more
prone to demonstrating responses which outlie from the expected and thus the
patients’ measurements may not be consistent with the population average and
change over time as they grow accustomed to the new treatments. Incorrectly
modelling such scenarios limits our ability to fully decipher the relationship
between how patients’ responses change over time and the impact this has on
their risk of an adverse event.

This work explores two examples, one on dialysis, and another on liver
cirrhosis. The first example analyses the haemoglobin levels of Northern
Ireland (NI) renal patients after commencing haemodialysis treatment. The
second example, utilising the Mayo Clinic Primary Biliary Cirrhosis data-set,
analyses data obtained from patients with primary biliary cirrhosis who
participated in a randomized placebo controlled trial of the drug
D-penicillamine. In both case studies, it is reasonable to assume that
patients may take time to adjust to the new treatments and thus would be more
prone to outlying responses.

A more detailed discussion on the two motivational examples is given in
Section 2. The rest of this paper is organised as follows. In Section 3, the
methodology for the various joint modelling approaches investigated in the
paper is presented, including the introduction of a robust joint model with
time-varying degrees-of-freedom. Section 4 presents a simulation study
comparing the approaches assuming an underlying structure of time-varying
degrees-of-freedom. Analysis of the two motivating examples is presented in
Section 5 followed by discussions in Section 6.

##  2 Motivating Examples

###  2.1 NI Renal data

Recent renal research has shown the potential for patients’ haemoglobin (Hb)
levels as an emerging biomarker for chronic kidney disease 15. It is known
within the renal research that in the initial stages after commencing
haemodialysis, patients’ haemoglobin levels tend to fluctuate to a large
degree and is deterimental to their survival. This research aims to capture
and model these initial stages of variability where there is a greater
potential for patients to either have outlying observations from their own
trends over time or be an outlier themselves from the population.

This work will analyse data that has been collected from seven different renal
centres by the NI Renal Information Service. The NI renal data-set contains
longitudinal information on 1,339 patients who have undergone haemodialysis
treatment between April 2002 and July 2011 4, 10, 16. In total, 27,064 Hb
measurements were recorded, with a median of 14 observations per individual.
Other health related measurements, including creatinine and urea levels, were
additionally monitored. Individuals were observed for a minimum of 2 days, and
a maximum of approximately 9 years throughout the study follow-up period. 44%
of the patients died during the study period with the remaining patients being
right-censored.

Figures 1 and 2 show the tail behaviour of the standardised conditional
residuals and random-effects, respectively, based on a joint model fitted for
the NI renal data with Gaussian assumptions for both random-effects and error
terms. Degrees-of-freedom estimates for Figure 1 were obtained by fitting a
univariate $t$ distribution such that $t(\mu,\sigma,\nu)$, where $\mu$ is the
location parameter, $\sigma$ the scale, $\nu$ the degrees-of-freedom. Time (in
years) was divided into bins based on 20% percentiles of the follow-up time.
As illustrated by Figures 1 and 2, there are departures from Normality for
both random-effects and error terms. Figure 1 further indicates that the
degrees-of-freedom parameter of the $t$-distribution for the residuals varies
over time, an anticipated result based on the known impact on patients’ Hb
levels initially as they adjust to the commencement of haemodialysis
treatment. Such changes over time will be accounted for through the time-
varying degrees-of-freedom formulation introduced in Section 3.4.

###  2.2 Primary biliary cirrhosis (PBC) data

A second motivating example arises from data collected from 312 patients with
primary biliary cirrhosis (PBC), a rare autoimmune disease which results in
cirrhosis of the liver and fatality. Patients suffering from PBC produce
increased levels of alkaline in the blood. The PBC data-set was collected from
the Mayo Clinic trial conducted between 1974 and 1984, which aimed to assess
the effectiveness of the D-penicillamine drug 17. Of the 312 patients, 158
were randomly assigned to the drug, and 154 to placebo. Median baseline age
was 48.8 (minimum = 26.28, maximum = 78.44). In total, 1,885 repeated alkaline
observations were collected over the 10 year trial period, with a median of 5
(minimum = 1, maximum = 16) observations per individual. Individuals were
observed for a minimum of 0.1 years, and a maximum of 14.3 years. 140 (45%)
patients died during the study-period.

The quantile-quantile plot for the standardised conditional residuals under
the joint model with Gaussian assumption (left panel of Figure 3) indicates
heavier tails than Gaussian. Estimated degrees-of-freedom from univariate t
distributions fitted to residuals that fall into bins based on 20% percentiles
of follow-up time indicate time-varying degrees-of-freedom structure up to
some extent. Figure 4 indicates that normality may be a reasonable assumption
for the random effects.

##  3 Approaches for Robust Joint Modelling

###  3.1 Notation

Before introducing the models, we present the general notation used thoughout
the paper. Let $Y_{ij}$ denote the $j^{\text{th}}$ $(j=1,\ldots,m_{i})$
repeated measurement belonging to subject $i$ $(i=1,\ldots,n)$ collected at
time $t_{ij}$, $t_{i}=\\{t_{i1},\ldots,t_{im_{i}}\\}$ the set of the follow-up
times at which $Y_{ij}$’s are collected, $a_{i}=\\{a_{i1}\ldots a_{1l}\\}$
baseline covariate information, and $S_{i}$ survival time. Here, $S_{i}$ is
subject to right-censoring. Hence, we introduce an additional random variable,
$E_{i}$, defined by $E_{i}=I(S_{i}^{*}\leq C_{i})$, where $I(\cdot)$ being
indicator function, $S_{i}^{*}$ true survival time for subject $i$, $C_{i}$
censoring time for the subject. $C_{i}$ is defined as $C_{i}=\min(C,D_{i})$,
with $C$ is the study end-time and $D_{i}$ is the drop-out time for subject
$i$.

###  3.2 Joint modelling of longitudinal and survival data

The framework for the so-called joint models for longitudinal and survival
data under the shared-parameter paradigm is given by

$\displaystyle Y_{ij}$ $\displaystyle=Y_{i}^{*}(t_{ij})+Z_{ij},$
$\displaystyle={\boldsymbol{x}}_{ij}^{\top}{\boldsymbol{\alpha}}+{\boldsymbol{d}}_{ij}^{\top}{\boldsymbol{B}}_{i}+Z_{ij},$
(1) $\displaystyle h_{i}(t)$
$\displaystyle=h_{0}(t)\exp\left({\boldsymbol{c}}_{i}^{\top}{\boldsymbol{\omega}}+f(\mathcal{Y}_{i}^{*}(t);{\boldsymbol{\eta}})\right).$
(2)

The framework allows for biomarker values (the longitudinal data) to be
measured with error, i.e. the observed data, $Y_{ij}$, is composed of
underlying continuous-time signal at time $t_{ij}$, $Y_{i}^{*}(t_{ij})$, and
noise, $Z_{ij}$. The signal is de-composed into fixed-effects,
${\boldsymbol{x}}_{ij}^{\top}{\boldsymbol{\alpha}}$, and random-effects,
${\boldsymbol{d}}_{ij}^{\top}{\boldsymbol{B}}_{i}$. Here,
${\boldsymbol{x}}_{ij}$ is a $p\times 1$ matrix structured by $a_{i}$ and
$t_{i}$. ${\boldsymbol{\alpha}}$ is a $p\times 1$ matrix of regression
coefficients as in multiple linear regression. ${\boldsymbol{d}}_{ij}$ is a
$q\times 1$ matrix, that is typically a subset of ${\boldsymbol{x}}_{ij}$.
${\boldsymbol{B}}_{i}$ are subject-specific coefficients that take into
account heterogeneity between subjects. $h_{i}(t)$ is the hazard of survival
event for subject $i$ at time $t$. $h_{0}(t)$ is the baseline hazard that can
be specified using hazard functions of parametric distributions, e.g. Weibull;
left un-specified as in Cox 18; assumed to be piece-wise constant; or be
expressed in terms of splines, e.g. natural cubic or B-splines.
${\boldsymbol{c}}_{i}$ is a $g\times 1$ matrix with elements from $a_{i}$,
${\boldsymbol{\omega}}$ $g\times 1$ vector of regression coefficients.
$f(\mathcal{Y}_{i}^{*}(t);{\boldsymbol{\eta}})$ is the term for taking into
account the association between hazard of survival event and features of
biomarker process, with $f(\cdot)$ being a known function. A popular choice is
to use the current value parametrisation such that
$f(\mathcal{Y}_{i}^{*}(t);{\boldsymbol{\eta}})=\eta Y_{i}^{*}(t)$. For other
choices, see Hickey et al. 19.

The classical model assumptions for the random-effects and error components
are zero-mean Gaussian distributions such that

$\displaystyle\boldsymbol{\mathrm{B}}_{i}$
$\displaystyle\sim\mathcal{N}\left(\boldsymbol{\mathrm{0}},\boldsymbol{\mathrm{\Sigma}}\right),$
$\displaystyle Z_{ij}$
$\displaystyle\sim\mathcal{N}\left(0,\sigma^{2}\right),$
$\displaystyle\boldsymbol{\mathrm{B}}_{i}\perp Z_{ij},\ \ $ $\displaystyle
Z_{ij}\perp Z_{ij^{\prime}}\ \mbox{for}\ j\neq j^{\prime},$

where $\boldsymbol{\mathrm{\Sigma}}$ is the covariance matrix of the random-
effects, and $\sigma$ is the standard deviation of the measurement error. In
what follows, we extend the assumptions beyond the Gaussian.

###  3.3 Robust joint modelling

The robustness of the joint model will be determined by the distributions of
the mixing variables, $\boldsymbol{\mathrm{B}}_{i}$ and $Z_{ij}$. The tail of
the density for ${\boldsymbol{\mathrm{B}}}_{i}$ determines robustness against
$b$-outliers and the tail of the density for $Z_{ij}$ determines the
robustness against $e$-outliers. We consider symmetric robust distributions
for both the random-effects and error components through normal-variance
mixtures such that

$\displaystyle\boldsymbol{\mathrm{B}}_{i}$
$\displaystyle=\boldsymbol{\mathrm{\Sigma}}^{1/2}\sqrt{V^{B}_{i}}\boldsymbol{\mathrm{B}}^{*}_{i},$
$\displaystyle Z_{ij}$ $\displaystyle=\sigma\sqrt{V^{Z}_{ij}}Z^{*}_{ij},$

where
$\boldsymbol{\mathrm{B}}^{*}_{i}\sim\mathcal{N}(\boldsymbol{\mathrm{0}},\boldsymbol{\mathrm{I}}_{q\times
q})$ and $Z^{*}_{ij}\sim\mathcal{N}(0,1)$. This formulation is flexible and
includes widely used distributions as special cases 20. The tail behaviour is
determined by the distribution of $V_{i}^{B}$ and $V_{ij}^{Z}$. The special
case $V_{i}^{B}=1$ and $V_{ij}^{Z}=1$ recovers the Gaussian joint model (see
Section 3.2). In this study, we specifically consider inverse Gamma
distribution $\mathcal{IG}$ with equal shape and scale parameters for the $V$
terms which results a t-distribution.

Ideally, one can use the posterior distribution of $V_{i}^{B}$ or $V_{ij}^{Z}$
to detect outliers. If the posterior distribution of $V_{i}^{B}$
($V_{ij}^{Z}$) is concentrated at large values this should indicate a
$b$-outlier for subject $i$ ($e$-outlier for observation $j$ of patient $i$).
However, in much of the robust modelling literature dependence between
$V_{i}^{B}$ and $V_{ij}^{Z}$, e.g. see 21, and $V_{ij}^{Z}$ and
$V_{ij^{\prime}}^{Z}$, e.g. Baghfalaki, Ganjali and Berridge 9, have been
introduced (mainly to simplify the inferential procedure). These dependencies,
hovewer, make the aforementioned interpretation on outlier detection
impossible.

In what follows, we will set three robust joint model formulations that are
available in the literature and discuss implications.

####  3.3.1 Approach 1

McCrink 10 considered

$\displaystyle V^{B}_{i}$ $\displaystyle=V_{i},$ $\displaystyle V^{Z}_{ij}$
$\displaystyle=V_{i},$ $\displaystyle V_{i}\sim$
$\displaystyle\mathcal{IG}(\gamma/2,\gamma/2).$

The mixing variable is equal for both random-effects and error components.
This implies that $b$\- and $e$-outliers must be present simultaneously for a
subject, which would be a strong assumption. Note that under this assumption
${\boldsymbol{\mathrm{B}}}_{i}$ and
${\boldsymbol{\mathrm{Z}}}_{i}=[Z_{i1}\ldots Z_{im_{i}}]^{\top}$ will be
jointly multivariate $t$ with a single degrees-of-freedom parameter, and the
following properties do not hold: $\boldsymbol{\mathrm{B}}_{i}\perp Z_{ij}$,
and $Z_{ij}\perp Z_{ij^{\prime}}$ for $j\neq j^{\prime}$.

####  3.3.2 Approach 2

Baghfalaki, Ganjali and Berridge 9 considered a robust joint model by setting

$\displaystyle V^{B}_{i}$
$\displaystyle=V^{B}_{i}\sim\mathcal{IG}(\phi/2,\phi/2),$ $\displaystyle
V^{Z}_{ij}$ $\displaystyle=V^{Z}_{i}\sim\mathcal{IG}(\delta/2,\delta/2).$

Under this approach, $\boldsymbol{\mathrm{B}}_{i}$ and
${\boldsymbol{\mathrm{Z}}}_{i}$ will be multivariate $t$, with separate
degrees-of-freedom parameters; the property of
$\boldsymbol{\mathrm{B}}_{i}\perp Z_{ij}$ holds, whereas $Z_{ij}\perp
Z_{ij^{\prime}}$ for $j\neq j^{\prime}$ does not. As $Z_{ij}$ and
$Z_{ij^{\prime}}$ for $j\neq j^{\prime}$ share a common $V^{Z}_{i}$, this
approach can be seen as a random-effects approach on the variance of $Z$. This
approach forces all the elements of $\\{Z_{ij}:j=1,\ldots,m_{i}\\}$ to be
outliers.

####  3.3.3 Approach 3

Asar, Fournier and Dantan 11 considered the setting of

$\displaystyle V^{B}_{i}$
$\displaystyle=V^{B}_{i}\sim\mathcal{IG}(\phi/2,\phi/2),$ $\displaystyle
V^{Z}_{ij}$ $\displaystyle=V^{Z}_{ij}\sim\mathcal{IG}(\delta/2,\delta/2),$

for joint modelling. Under this approach the dependence between $Z_{ij}$ and
$Z_{ij^{\prime}}$ that was present in Approach 2 has been removed. This is the
most natural approach compared to the previous two, and will be used as a base
for our time varying degrees-of-freedom formulation.

###  3.4 Appoach 4: Time-varying degrees-of-freedom formulation

By building on Approach 3, we consider time-varying degrees-of-freedom
parameter, $\delta(t)$, for $Z_{ij}$. In practice, we need to discretise time,
i.e. $\delta(t)=\delta(t_{ij})$. This approach considers

$\displaystyle V^{Z}_{ij}$
$\displaystyle\sim\mathcal{IG}(\delta_{ij}/2,\delta_{ij}/2),$
$\displaystyle\delta_{ij}$
$\displaystyle=\exp(\delta_{0}+{\boldsymbol{a}}_{ij}^{\top}{\boldsymbol{\beta}}).$

Here, ${\boldsymbol{a}}_{ij}$ is a $s\times 1$ to be specified by natural
cubic splines or B-splines, ${\boldsymbol{\beta}}$ associated coefficients.
Note that one obtains time invariant degrees-of-freedom as
$\delta_{ij}=\delta=\exp(\delta_{0})$, when
${\boldsymbol{\beta}}={\boldsymbol{0}}$.

###  3.5 Priors

We only give details about the prior distributions that are assigned to the
parameters based on Approach 4, since Approaches 1-3 are just special cases.
We set weakly-informative prior distributions for the parameters. Elements of
$\alpha$ are assigned zero-mean Cauchy prior with scale of 5,
$\mathcal{C}(0,5)$, only $\alpha_{0}$ was given $\mathcal{C}(0,20)$.
${\boldsymbol{\mathrm{\Sigma}}}$ is re-written as
${\boldsymbol{\mathrm{R}}}{\boldsymbol{\mathrm{\Omega}}}{\boldsymbol{\mathrm{R}}}$,
with ${\boldsymbol{\mathrm{R}}}$ being a diagonal matrix of scale parameters
of $B_{h}$ $(h=1,\ldots,q)$, and ${\boldsymbol{\mathrm{\Omega}}}$ a
correlation matrix 22. Elements of ${\boldsymbol{\mathrm{R}}}$ are given half-
Cauchy, $\mathcal{C}_{+}(0,5)$, whereas elements of
${\boldsymbol{\mathrm{\Omega}}}$ are given LKJ with the parameter of 2,
$\mathcal{LKJ}(2)$. Degrees-of-freedom parameter for
${\boldsymbol{\mathrm{B}}}_{i}$, $\phi$, and $\delta_{0}$ of $Z_{ij}$ are
given uniform priors, between 2 and 100. Elements of
${\boldsymbol{\mathrm{\beta}}}$ are given $\mathcal{C}(0,5)$. $\sigma$ is
given $\mathcal{C}_{+}(0,5)$. Log-transformed elements of $h_{0}(t)$ and
elements of ${\boldsymbol{\mathrm{\omega}}}$ and
${\boldsymbol{\mathrm{\eta}}}$ are given $\mathcal{C}(0,5)$.

###  3.6 Inference

For inference, we use Markov Chain Monte Carlo methods to sample from the
joint posterior densities of the parameters and latent variables given data.
We specifically consider the Hamiltonian Monte Carlo (HMC) 23 as implemented
in the No-U-Turn Sampler (NUTS) 24. We do not present the details of neither
the likelihood function, nor HMC and NUTS, since, whereas the former is quite
straightforward (e.g. see 11), the second can be followed from the cited
references.

Bespoke R 25 codes to fit the joint models under Approaches 1-4 are available
from the robjm package (https://github.com/ozgurasarstat/robjm) that
internally uses the HMC sampling engine Stan 26 through the R package Rstan
27.

##  4 Simulation study

A simulation study was conducted to investigate the effects of time-varying
degrees-of-freedom in the estimation of parameters for the joint modelling
approaches discussed in Section 3. A sample size of $n=250$ individuals was
considered with 200 data-sets being simulated under the assumption of time-
varying degrees-of-freedom for the longitudinal residuals. Note that we
consider 200 as the number of replications mainly because of the computational
cost.

Data were generated under Approach 4. The assumed underlying structure for the
time-varying degrees-of-freedom is illustrated in Figure 5. The longitudinal
sub-model was given by:

$\displaystyle Y_{ij}$ $\displaystyle=$ $\displaystyle
Y^{*}_{i}(t_{ij})+Z_{ij},$ (3) $\displaystyle=$
$\displaystyle\alpha_{1}+t_{ij}\alpha_{2}+B_{1i}+t_{ij}B_{2i}+Z_{ij},$

where each individual has an average of 20 observations between time points 0
and 5. A random-intercept ($\alpha_{1}+B_{1i}$) and random-slope
($\alpha_{2}+B_{2i}$) model was assumed to replicate what is most commonly
used in the joint modelling literature. Survival data was generated from the
following model:

$\displaystyle h_{i}(t)$
$\displaystyle=h_{0}(t)\exp\left(X_{i}\omega+Y_{i}^{*}(t)\eta\right),$ (4)

where $X\sim\textrm{Bernoulli}(0.5)$, $h_{0}(t)$ specified by a Weibull
baseline hazard, i.e. $h_{0}(t)=\lambda\nu t^{\nu-1}$. There is a final
truncation time of 5 after which non-informative right-censoring occurs. The
true values for the unknown parameters are given in Table 1.

Each robust model and the model with Gaussian assumptions for ${\bf B}_{i}$
and $Z_{ij}$, called the standard joint model, were estimated utilising the
robjm package with four chains each of length 2,000 with the first 1000
iterations considered as warm-up. The averages of posterior means, width of
the corresponding 95% credibility intervals (CI) and coverage (Cov.) of 200
replications for each model are collated in Table 1.

For all the models, including the standard joint model, fixed effects for both
the longitudinal and survival models and baseline hazard parameters
demonstrate similar biases and coverages. In terms of averages of the
credibility intervals for $\alpha_{1}$ and $\alpha_{2}$, the standard joint
model produces larger values compared to the other models. Coverages for the
elements of ${\bf\Sigma}$ are lower than expected for the standard joint
model, and almost 0 for Approach 1, whereas they are at the expected level for
Approaches 2, 3 and 4. Approaches 1 and 2 produces biased and 0 coverage
results for $\sigma^{2}$. This highlights the need to remove the dependence
between $Z_{ij}$ and $Z_{ij^{\prime}}$. Approach 3 provides similar estimates
and levels of coverage compared to Approach 4. Therefore, whilst Approach 3
has the inability to fully capture the time-varying nature of the degrees-of-
freedom for the residuals, this appears to have limited impact on the
estimation of the other parameters.

##  5 Illustrative Data Examples

###  5.1 Renal Example

The first case-study considers NI haemodialysis patients’ data with a joint
model specified as:

$\displaystyle Y_{ij}$ $\displaystyle=$ $\displaystyle Y_{i}(t_{ij})+Z_{ij},$
(5) $\displaystyle=$
$\displaystyle\alpha_{1}+t_{ij}\alpha_{2}+\text{age}_{i}\alpha_{3}+\text{creatinine}_{ij}\alpha_{4}+\text{MCHC}_{ij}\alpha_{5}+\text{MCV}_{ij}\alpha_{6}+\text{urea}_{ij}\alpha_{7}$
$\displaystyle+\text{male}_{i}\alpha_{8}+\text{I(EPO}_{ij}\text{=Aranesp)}\alpha_{9}+\text{I(EPO}_{ij}\text{=Alfa)}\alpha_{10}+\text{I(EPO}_{ij}\text{=Other)}\alpha_{11}$
$\displaystyle+\text{Prior2007}_{i}\alpha_{12}+\text{IronHydroxide}_{ij}\alpha_{13}+\text{Venofer}_{ij}\alpha_{14}+B_{1i}+t_{ij}B_{2i}+Z_{ij}$
$\displaystyle h_{i}(t)$ $\displaystyle=$ $\displaystyle
h_{0}(t)\exp\left(\text{age}_{i}\omega+Y_{i}^{*}(t)\eta\right),$ (6)

where the response $Y$ is an individual’s Hb level, $h_{0}(t)$ is assumed to
be Weibull baseline hazard, MCHC represents mean corpuscular haemoglobin
concentration, MCV represents mean corpuscular volume and EPO represents
erythropoietin drug treatments.

The models were fitted using the robjm package with four chains each of length
8,000 with the first 2,000 iterations considered as warm-up. Posterior
summaries for each of the approaches are given within Table 2. Convergence was
confirmed using traceplots and R-hat statistics 28.

Similar to the findings of the simulation study in Section 4, in general,
results based on Approaches 3 and 4 are quite similar. There are differences
in the results of these two models and the others, some parameters were
overestimated whereas some were underestimated. For example, the association
parameter, $\eta$, appears to be affected by the choice of model. It is
underestimated by the standard joint model, which was overcome by Approach 1
up to some extent. Smoothed density plots for $\eta$ are displayed in the left
panel of Figure 6. Posterior summaries of degrees-of-freedom for Approach 3
and Approach 4 with 7 knots are displayed in the right panel of Figure 6.

###  5.2 PBC Example

The second case-study uses data from the PBC study. The model we fit has the
following form:

$\displaystyle Y_{ij}$ $\displaystyle=Y_{i}^{*}(t_{ij})+Z_{ij}$
$\displaystyle=\alpha_{1}+\text{age}_{i}\alpha_{2}+t_{ij}\alpha_{3}+\text{I(drug${}_{i}$
= D-penicillamine)}\alpha_{4}+B_{i1}+t_{ij}B_{i2}+Z_{ij},$ (7) $\displaystyle
h_{i}(t)$
$\displaystyle=h_{0}(t)\exp\left(\text{age}_{i}\omega_{1}+\text{I(drug${}_{i}$
= D-penicillamine)}\omega_{2}+Y_{i}^{*}(t)\eta\right),$ (8)

where $Y=\log(\text{alkaline})$, age is baseline age, $t$ follow-up time,
$h_{0}(t)$ Weibull baseline hazard function.

For Approaches 2–4, the results for $\phi$ indicated that Normal assumption
might be reasonable for ${\boldsymbol{\mathrm{B}}}_{i}$. Therefore, only $Z$
term was assumed to be $t$-distributed. For Approach 4, we considered number
of knots ranging from 1 to 5. The knots were put into the empirical quantiles
of the time variable. Four chains, each length of 6,000, were started from
random initials. Half of each chain was considered as warm-up. Convergence was
checked through traceplots and R-hat statistics 28. Posterior summaries,
specifically the 2.5%, 50% and 97.5% percentiles, are displayed in Table 3.

The association parameter, $\eta$, seems to be affected by the model choice.
Regarding model comparisons in terms of parameter estimates, similar
observations to those made for the NI renal application were observed for this
case-study. The smoothed densities of the MCMC samples of $\eta$ and the
degrees-of-freedom estimates for Approach 3 and Approach 4 with 1 knot are
displayed in the left and right panels of Figure 7, respectively.

R codes to run the models for the PBC data could be reach at
https://github.com/ozgurasarstat/robjm/blob/master/tests/testthat/test_pbc_analysis.R.

##  6 Discussion

This work introduces a new flexible approach to model longitudinal outliers in
a joint modelling framework, where the degrees-of-freedom parameter for the
residuals is assumed to vary over time. This scenario replicates the common
situation when patients take time to adjust to new treatments, resulting in
more outlying measurements (e-outliers) at different times across the period
of observation. In addition, this paper contrasts the new time-varying
approach with three alternative time-invariant formulations of robust joint
models currently found in the literature and the standard joint model with
Gaussian assumptions for the random terms. Each of the approaches presented
can be fitted in the accompanying robjm software package.

Throughout the results presented within the simulation study and the real-life
data analysis, the need to properly account for longitudinal outliers is
evident, a practice not widely adopted in the current joint modelling
literature. This paper highlights the need to investigate the presence of
outliers and provides a flexible approach to downweighing their negative
impact, contrasting the effect of various assumptions made by the different
approaches to robust joint modelling. It is evident that bias and inefficiency
results from ignorance of longitudinal outliers and the assumption of
normality for the random terms. In addition, this work highlights the need to
remove the restrictive assumptions of Approaches 1 and 2.

Unexpectedly, despite a true underlying assumption of time-varying degrees-of-
freedom within the simulation study, Approach 3 was found to be in good
agreement with the time-varying model. Within the simulation study and in the
two real-life examples, the degrees-of-freedom were assumed or estimated to be
consistently low and varying over a limited range throughout the period of
observation. Under such scenarios, Approach 3 appears to provide a time-
invariant average of the time-varying degrees-of-freedom for the residuals
over the period of observation, as demonstrated within Figure 6 for the renal
analysis. Given the limited change in estimated parameters between the two
approaches, it is not therefore unexpected to find high levels of agreement
when comparing the posterior summaries of the predictions for the longitudinal
random effects, as illustrated in Figure 8 from the PBC analysis.

This work highlights the robustness provided by Approach 3 in the estimation
of the degrees-of-freedom, providing similar estimates as the time-varying
approach inspite of the presence of time-varying degrees-of-freedom. It
signifies that in cases where the degrees-of-freedom parameter does not vary
widely, Approach 3 is recommended when longitudinal outliers are present.

We did not consider model selection methods to select the number of knots,
mainly because the number of knots does not seem to change the results
considerably. In the time-varying degrees-of-freedom model we assumed the
scale parameter, $\sigma$ being time-constant. This parameter might also be
assumed to be time-varying.

## Funding

This work was supported by the Engineering and Physical Sciences Research
Council [Reference: EP/P026028/1]

### Acknowledgment

Dr. Laura Boyle’s (The University of Adelaide) helps are gratefully
acknowledged.

### Conflict of interest

The authors declare no potential conflict of interests.

## References

  * Tsiatis and Davidian 2004 Tsiatis A.A. and Davidian M. (2004). Joint modeling of longitudinal and time-to-event data: An overview. _Statistica Sinica_ 14(3), 809-834.
  * Diggle, Sousa and Chetwynd 2008 Diggle P.J., Sousa I. and Chetwynd A.G. (2008). Joint modelling of repeated measurements and time-to-event outcomes: The fourth Armitage lecture. _Statistics in Medicine_ 27(16), 2981-2998.
  * Rizopoulos 2012 Rizopoulos D. (2012). Joint Models for Longitudinal and Time-To-Event Data: With Applications in R. Chapman & Hall/CRC Press.
  * McCrink, Marshall and Cairns 2013 McCrink L.M., Marshall A.H. and Cairns K.J. (2013). Advances in Joint Modelling: A review of recent developments in software and techniques using the modelling of survival of end stage renal disease patients as a working example. _International Statistical Review_ 81, 249-269.
  * Asar et al. 2015 Asar Ö., Ritchie J., Kalra P. A. and Diggle P. J. (2015). Joint modelling of repeated measurement and time-to-event data: an introductory tutorial. _International Journal of Epidemiology_ 44(1), 334-344.
  * Elashoff 2017 Elashoff R.M. (2017). Joint modeling of longitudinal and time-to-event data. CRC Press.
  * Li, Elashoff and Li 2009 Li N., Elashoff R.M. and Li G. (2009). Robust joint modeling of longitudinal measurements and competing risks failure time data. _Biometrical Journal_ 51(1), 19-30.
  * Huang, Li and Elashoff 2010 Huang X., Li G. and Elashoff R. (2010). A joint model of longitudinal and competing risks survival data with heterogeneous random effects and outlying longitudinal measurements. _Statistics and Its Interface_ 3(2), 185-195.
  * Baghfalaki, Ganjali and Berridge 2013 Baghfalaki T., Ganjali M. and Berridge D. (2013). Robust joint modeling of longitudinal measurements and time to event data using normal/independent distributions: A Bayesian approach. _Biometrical Journal_ 34(2), 187-220.
  * McCrink 2014 McCrink L.M. (2014). Outlier effects on robust joint modelling of longitudinal and survival data. _PhD Thesis, Queen’s University Belfast_.
  * Asar, Fournier and Dantan 2019 Asar Ö., Fournier M.-C. and Dantan E. (2019). Dyanmic predictions of kidney graft survival in the presence of longitudinal outliers. _Under review_. Avaiable at https://arxiv.org/abs/1905.00816.
  * Taylor et al. 2013 Taylor J.M.G., Park Y., Ankerst D.P., Proust-Lima C., Williams S., Kestin L., Bae K., Pickles T., and Sandler H. (2013). Real-time individual predictions of prostate cancer recurrence using joint models. _Biometrics_ 69(1), 206-213.
  * Song, Peng and Tu 2012 Song H., Peng Y., and Tu D. (2012). A new approach for joint modelling of longitudinal measurements and survival times with a cure fraction. _Canadian Journal of Statistics_ 40(2), 207-224.
  * Baghfalaki, Ganjali and Hashemi 2014 Baghfalaki T., Ganjali M. and Hashemi R. (2014). Bayesian Joint Modeling of Longitudinal Measurements and Time-to-Event Data Using Robust Distributions. _Journal of Biopharmaceutical Statistics_ 24(4), 834-855.
  * Chen et al. 2015 Chen Y., Qin M., Zheng J., Yan H., Li M., Cui Y., Zhang R., Zhao W. and Guo Y. (2015). Hemoglobin discriminates stages of chronic kidney disease in elderly patients. _Experimental and Therapeutic Medicine_ 10, 567-571.
  * Donnelly et al. 2017 Donnelly, C., McFetridge, L. M., Marshall, A. H. and Mitchell, H. J. (2017). A two-stage approach to the joint analysis of longitudinal and survival data utilising the Coxian phase-type distribution. _Statistical Methods in Medical Research_ 27(12), 3577-3594.
  * Murtaugh et al. 1994 Murtaugh P., Dickson E., Van Dam G., Malincho M., Grambsch P., Langworthy A., and Gips C. (1994) Primary biliary cirrhosis: prediction of short-term survival based on repeated patient visits. _Hepatology_ 20, 126-134.
  * Cox 1972 Cox D. R. (1972). Regression models and life-tables (with discussion). _Journal of the Royal Statistical Society, Series B (Methodological)_ 34(2), 187-220.
  * Hickey et al. 2016 Hickey G. L., Philipson P., Jorgensen A. and Kolamunnage-Dona R. (2016). Joint modelling of time-to-event and multivariate longitudinal outcomes: recent developments and issues. _BMC Medical Research Methodology_ 16, 117\.
  * Asar et al. 2018 Asar Ö., Bolin D., Diggle P. J. and Wallin J. (2018). Linear mixed-effects models for non-gaussian repeated measurement data. _Under review_. Avaiable at arxiv.org/abs/1804.02592.
  * Pinheiro, Liu and Wu 2001 Pinheiro J. C., Liu C., and Wu Y. N. (2001). Efficient algorithms for robust estimation in linear mixed-effects models using the multivariate t distribution. Journal of Computational and Graphical Statistics 10, 249–276.
  * Stan Development Team 2019 Stan Development Team (2019). Stan User’s Guide, Version 2.19. URL mc-stan.org/docs/2_19/stan-users-guide-2_19.pdf.
  * Neal 2011 Neal R. (2011). MCMC Using Hamiltonian Dynamics in Handbook of Markov Chain Monte Carlo, pages 113-162. Eds. Brooks S., Gelman A., Jones G. L. and Meng, X. L. Chapman & Hall/CRC Press.
  * Hoffman and Gelman 2014 Hoffman, M. D. and Gelman, A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo _Journal of Machine Learning Research_ 15, 1593-1623.
  * R Core Team 2018 R Core Team(2018). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
  * Carpenter et al. 2017 Carpenter B., Gelman A., Hoffman M. D., Lee D., Goodrich B., Betancourt M., Brubaker M., Guo J., Li P. and Riddell A. (2017). Stan: A Probabilistic Programming Language. _Journal of Statistical Software_ 1, 1-32.
  * Stan Development Team 2018 Stan Development Team (2018). RStan: the R interface to Stan. R package version 2.17.3. URL http://mc-stan.org/.
  * Brooks and Gelman 1997 Brooks S. P. and Gelman A. (1997). General methods for monitoring convergence of iterative simulations. _Journal of Computational and Graphical Statistics_ , 7, 434–455.

Figure 1: Plots to inspect tail behavior of the standardised conditional
residuals obtained from the joint model with Gaussian assumptions for the
renal dataset. Details of calculations are provided within Section 2.1. Left:
Quantile-quantile plot against standard Normal distribution; Right: Degrees-
of-freedom parameter estimates. Figure 2: Plots to inspect the distributions
of the random-intercept (left) and random-slope (right) terms against standard
Normal based on the renal data. Figure 3: Plots to inspect the standardised
conditional residuals obtained from the joint model with Gaussian assumptions
for the PBC data-set. Left: Quantile-quantile plot against standard Normal
distribution, Right: degrees-of-freedom parameter estimates. Degrees-of-
freedom are obtained similar to Figure 1. Figure 4: Plots to inspect the
distributions of the random-intercept (left) and random-slope (right) terms
against standard Normal. Figure 5: Underlying degrees-of-freedom assumption
for simulation study. Figure 6: Left panel: Posterior densities of $\eta$
under different models. Right panel: Medians of degrees-of-freedom parameters
(mid-lines) and 2.5% and 97.5% percentiles of the MCMC samples under Approach
3 (solid) and time-varying model with 7 knots (dashed) based on NI renal data.
Figure 7: Left panel: smoothed MCMC samples of $\eta$ under different models.
Right panel: Medians of degrees-of-freedom parameters (mid-lines) and 2.5% and
97.5% percentiles of the MCMC samples under model 3 and tv model with 1 knot.
Figure 8: Posterior summaries (2.5%, 50% and 97.5% percentiles of the MCMC
samples) for Approach 3 against the time-varying model with 5 knots from the
PBC analysis for the predictions of the longitudinal random intercept (top
row) and slope (bottom row).

Table 1: Table showing the simulation results based on data assuming time-varying degrees-of-freedom |  |  |  | Standard Joint | Approach 1 | Approach 2 | Approach 3 | Approach 4
---|---|---|---|---|---|---|---|---
|  |  | True | Mean | CI | Cov. | Mean | CI | Cov. | Mean | CI | Cov. | Mean | CI | Cov. | Mean | CI | Cov.
|  | $\alpha_{1}$ | 1.00 | 0.992 | 0.329 | 0.930 | 0.986 | 0.291 | 0.935 | 0.994 | 0.247 | 0.920 | 0.995 | 0.247 | 0.910 | 0.995 | 0.246 | 0.920
|  | $\alpha_{2}$ | 0.40 | 0.391 | 0.233 | 0.940 | 0.389 | 0.209 | 0.950 | 0.396 | 0.175 | 0.935 | 0.396 | 0.176 | 0.930 | 0.396 | 0.175 | 0.930
|  | ${\boldsymbol{\Sigma}}(1,1)^{*}$ |  0.60 | 1.703 | 0.638 | 0.495 | 1.234 | 0.523 | 0.005 | 0.622 | 0.377 | 0.955 | 0.619 | 0.373 | 0.970 | 0.616 | 0.372 | 0.955
|  | ${\boldsymbol{\Sigma}}(2,1)^{*}$ |  0.25 | 0.673 | 0.363 | 0.590 | 0.497 | 0.287 | 0.070 | 0.251 | 0.184 | 0.960 | 0.251 | 0.183 | 0.955 | 0.252 | 0.183 | 0.955
|  | ${\boldsymbol{\Sigma}}(2,2)^{*}$ |  0.30 | 0.823 | 0.328 | 0.490 | 0.608 | 0.265 | 0.000 | 0.311 | 0.188 | 0.960 | 0.311 | 0.188 | 0.965 | 0.310 | 0.187 | 0.965
|  | ${\gamma}^{**}$ | — | — | — | — | 20.719 | — | — | — | — | — | — | — | — | — | — | —
|  | $\phi$ | 3.00 | — | — | — | — | — | — | 3.169 | 2.285 | 0.930 | 3.163 | 2.286 | 0.940 | 3.159 | 2.270 | 0.940
|  | ${\sigma^{2}}^{***}$ | 0.25 | 0.351 | — | — | 0.322 | 0.044 | 0.000 | 0.325 | 0.043 | 0.000 | 0.250 | 0.042 | 0.965 | 0.259 | 0.044 | 0.825
|  | ${\delta}^{***}$ | 5.00 | — | — | — | — | — | — | 30.350 | — | — | 7.149 | — | — | — | — | —
|  | $\delta_{0}$ | 5.00 | — | — | — | — | — | — | — | — | — | — | — | — | 5.801 | 6.160 | 0.950
|  | $\beta_{1}$ | 0.50 | — | — | — | — | — | — | — | — | — | — | — | — | 0.536 | 2.285 | 0.960
|  | $\beta_{2}$ | 0.50 | — | — | — | — | — | — | — | — | — | — | — | — | 0.668 | 2.648 | 0.945
|  | $\beta_{3}$ | 0.50 | — | — | — | — | — | — | — | — | — | — | — | — | 0.612 | 2.508 | 0.935
|  | $\beta_{4}$ | 0.50 | — | — | — | — | — | — | — | — | — | — | — | — | 0.614 | 2.41 | 0.935
|  | $\beta_{5}$ | 0.50 | — | — | — | — | — | — | — | — | — | — | — | — | 0.693 | 2.921 | 0.940
|  | $\beta_{6}$ | -0.25 | — | — | — | — | — | — | — | — | — | — | — | — | 0.152 | 2.400 | 0.895
|  | $\lambda$ | 0.04 | 0.040 | 0.036 | 0.950 | 0.040 | 0.036 | 0.955 | 0.040 | 0.036 | 0.950 | 0.040 | 0.036 | 0.945 | 0.040 | 0.036 | 0.950
|  | $\nu$ | 1.20 | 1.207 | 0.433 | 0.950 | 1.206 | 0.434 | 0.955 | 1.210 | 0.430 | 0.945 | 1.210 | 0.431 | 0.950 | 1.211 | 0.432 | 0.950
|  | $\omega$ | 0.50 | 0.494 | 0.709 | 0.945 | 0.494 | 0.709 | 0.945 | 0.494 | 0.708 | 0.940 | 0.494 | 0.709 | 0.935 | 0.494 | 0.713 | 0.950
|  | $\eta$ | 0.30 | 0.308 | 0.139 | 0.955 | 0.309 | 0.140 | 0.945 | 0.306 | 0.137 | 0.960 | 0.307 | 0.138 | 0.950 | 0.307 | 0.138 | 0.950
∗ For standard joint model, true values for these parameters were considered
as 1.8 $\left(0.6*\frac{3}{3-2}\right)$, 0.75 $\left(0.3*\frac{3}{3-2}\right)$
and 0.9 $\left(0.3*\frac{3}{3-2}\right)$.
∗∗ For this parameter, CI and Cov. were not calculated, since the truth is
unknown.
∗∗∗ For these parameters, CI and Cov. were not calculated, since the truths
are time-varying.

t  nor-nor tt-mod1 tt-mod2 tt-mod3 tv-2knot tv-3knot tv-4knot tv-5knot
tv-6knot tv-7knot $\alpha_{1}$ 10.138 10.239 10.261 10.246 10.245 10.248
10.247 10.251 10.249 10.250 (10.018, 10.259) (10.117, 10.354) (10.145, 10.376)
(10.133, 10.361) (10.130, 10.357) (10.132, 10.360) (10.134, 10.360) (10.135,
10.364) (10.132, 10.364) (10.136, 10.365) $\alpha_{2}$ -0.133 -0.138 -0.130
-0.127 -0.126 -0.127 -0.126 -0.126 -0.126 -0.126 (-0.163, -0.103) (-0.167,
-0.110) (-0.157, -0.103) (-0.154, -0.100) (-0.154, -0.099) (-0.154, -0.099)
(-0.154, -0.099) (-0.154, -0.099) (-0.154, -0.099) (-0.154, -0.099)
$\alpha_{3}$ 0.124 0.112 0.104 0.104 0.105 0.105 0.105 0.105 0.105 0.105
(0.089, 0.159) (0.077, 0.147) (0.071, 0.136) (0.073, 0.137) (0.073, 0.138)
(0.073, 0.137) (0.073, 0.137) (0.072, 0.137) (0.073, 0.138) (0.073, 0.138)
$\alpha_{4}$ 0.138 0.125 0.123 0.127 0.128 0.128 0.128 0.128 0.128 0.128
(0.121, 0.154) (0.109, 0.142) (0.107, 0.140) (0.111, 0.144) (0.111, 0.145)
(0.111, 0.144) (0.112, 0.144) (0.111, 0.144) (0.111, 0.144) (0.111, 0.145)
$\alpha_{5}$ 0.257 0.259 0.257 0.259 0.259 0.259 0.259 0.259 0.259 0.258
(0.239, 0.275) (0.242, 0.276) (0.240, 0.274) (0.241, 0.276) (0.241, 0.276)
(0.241, 0.276) (0.241, 0.276) (0.241, 0.276) (0.241, 0.276) (0.241, 0.276)
$\alpha_{6}$ 0.265 0.290 0.281 0.274 0.277 0.277 0.277 0.277 0.276 0.276
(0.222, 0.308) (0.245, 0.334) (0.238, 0.325) (0.231, 0.318) (0.234, 0.320)
(0.234, 0.321) (0.233, 0.319) (0.233, 0.320) (0.233, 0.320) (0.232, 0.320)
$\alpha_{7}$ 0.097 0.102 0.101 0.109 0.110 0.109 0.110 0.110 0.109 0.109
(0.075, 0.120) (0.080, 0.125) (0.079, 0.123) (0.087, 0.132) (0.087, 0.132)
(0.087, 0.132) (0.087, 0.132) (0.087, 0.132) (0.087, 0.131) (0.087, 0.132)
$\alpha_{8}$ -0.155 -0.142 -0.178 -0.187 -0.187 -0.190 -0.190 -0.191 -0.191
-0.191 (-0.255, -0.055) (-0.238, -0.045) (-0.270, -0.082) (-0.281, -0.095)
(-0.279, -0.095) (-0.280, -0.100) (-0.280, -0.097) (-0.282, -0.098) (-0.282,
-0.099) (-0.282, -0.100) $\alpha_{9}$ 0.207 0.187 0.176 0.195 0.193 0.193
0.192 0.192 0.191 0.191 (0.136, 0.278) (0.117, 0.257) (0.107, 0.243) (0.125,
0.261) (0.124, 0.261) (0.125, 0.263) (0.123, 0.261) (0.124, 0.258) (0.123,
0.261) (0.122, 0.258) $\alpha_{10}$ -0.037 -0.062 -0.064 -0.053 -0.051 -0.051
-0.050 -0.051 -0.052 -0.052 (-0.147, 0.074) (-0.171, 0.045) (-0.170, 0.041)
(-0.161, 0.054) (-0.160, 0.054) (-0.158, 0.055) (-0.157, 0.057) (-0.157,
0.054) (-0.160, 0.055) (-0.158, 0.057) $\alpha_{11}$ -0.100 -0.140 -0.152
-0.119 -0.131 -0.130 -0.130 -0.135 -0.139 -0.140 (-0.342, 0.151) (-0.381,
0.099) (-0.385, 0.089) (-0.346, 0.101) (-0.357, 0.092) (-0.357, 0.096)
(-0.356, 0.093) (-0.364, 0.088) (-0.366, 0.087) (-0.364, 0.089) $\alpha_{12}$
0.348 0.333 0.331 0.325 0.324 0.324 0.324 0.322 0.323 0.323 (0.246, 0.451)
(0.236, 0.433) (0.239, 0.426) (0.231, 0.419) (0.230, 0.417) (0.231, 0.418)
(0.228, 0.418) (0.206, 0.419) (0.228, 0.417) (0.230, 0.419) $\alpha_{13}$
0.329 0.292 0.282 0.310 0.310 0.310 0.311 0.311 0.311 0.311 (0.271, 0.388)
(0.236, 0.348) (0.227, 0.337) (0.255, 0.367) (0.254, 0.366) (0.254, 0.367)
(0.254, 0.367) (0.254, 0.367) (0.255, 0.367) (0.254, 0.368) $\alpha_{14}$
0.241 0.223 0.218 0.221 0.221 0.221 0.221 0.220 0.220 0.220 (0.201, 0.282)
(0.185, 0.261) (0.180, 0.255) (0.182, 0.260) (0.182, 0.259) (0.182, 0.260)
(0.182, 0.260) (0.179, 0.259) (0.181, 0.259) (0.182, 0.259) $\Sigma_{11}$
1.043 0.901 0.551 0.543 0.540 0.538 0.539 0.534 0.536 0.536 (0.940, 1.158)
(0.806, 1.005) (0.452, 0.664) (0.452, 0.648) (0.448, 0.644) (0.445, 0.640)
(0.446, 0.644) (0.442, 0.638) (0.447, 0.638) (0.444, 0.643) $\Sigma_{12}$
-0.256 -0.218 -0.136 -0.140 -0.138 -0.137 -0.137 -0.136 -0.137 -0.137 (-0.303,
-0.213) (-0.258, -0.182) (-0.172, -0.108) (-0.175, -0.110) (-0.173, -0.108)
(-0.171, -0.107) (-0.173, -0.108) (-0.171, -0.107) (-0.170, -0.108) (-0.172,
-0.107) $\Sigma_{22}$ 0.134 0.113 0.075 0.080 0.078 0.078 0.077 0.076 0.077
0.077 (0.111, 0.162) (0.094, 0.138) (0.059, 0.095) (0.063, 0.101) (0.061,
0.098) (0.061, 0.097) (0.061, 0.098) (0.060, 0.097) (0.061, 0.097) (0.061,
0.098) $\sigma^{2}$ 1.260 1.010 0.969 0.894 0.894 0.894 0.895 0.897 0.895
0.897 (1.238, 1.283) (0.972, 1.048) (0.931, 1.008) (0.866, 0.925) (0.865,
0.925) (0.866, 0.925) (0.866, 0.926) (0.868, 0.925) (0.866, 0.926) (0.866,
0.928) $\gamma$ — 8.685 — — — — — — — — (7.708, 9.802) $\phi$ — — 3.150 2.764
2.760 2.742 2.751 2.727 2.726 2.727 (2.504, 4.132) (2.280, 3.419) (2.273,
3.407) (2.263, 3.382) (2.283, 3.417) (2.244, 3.364) (2.261, 3.370) (2.256,
3.392) $\delta$ — — 7.789 6.847 — — — — — — (6.908, 8.804) (6.296, 7.561)
$\delta_{0}$ — — — — 5.532 5.071 4.671 4.894 4.554 4.355 (4.722, 6.600)
(4.167, 6.250) (3.782, 5.911) (3.874, 6.338) (3.553, 6.016) (3.347, 5.833)
$\beta_{1}$ — — — — 0.609 0.277 0.212 0.501 0.336 0.261 (0.272, 0.957) (0.041,
0.513) (-0.056, 0.489) (0.182, 0.833) (-0.032, 0.695) (-0.122, 0.652)
$\beta_{2}$ — — — — 0.089 0.704 0.589 0.044 0.546 0.694 (-0.346, 0.529)
(0.333, 1.097) (0.258, 0.921) (-0.348, 0.435) (0.106, 0.984) (0.234, 1.178)
$\beta_{3}$ — — — — -0.865 0.281 0.432 0.781 0.129 0.233 (-1.363, -0.328)
(-0.249, 0.830) (0.006, 0.858) (0.447, 1.115) (-0.263, 0.513) (-0.181, 0.667)
$\beta_{4}$ — — — — — -0.991 0.468 0.165 0.983 0.542 (-1.548, -0.379) (-0.165,
1.074) (-0.288, 0.619) (0.600, 1.396) (0.086, 0.984) $\beta_{5}$ — — — — — —
-0.812 0.198 -0.049 0.939 (-1.407, -0.162) (-0.483, 0.879) (-0.548, 0.459)
(0.526 ,1.340) $\beta_{6}$ — — — — — — — -0.526 0.399 -0.004 (-1.156, 0.204)
(-0.317, 1.132) (-0.516, 0.499) $\beta_{7}$ — — — — — — — — -0.358 0.524
(-1.052, 0.442) (-0.230, 1.263) $\beta_{8}$ — — — — — — — — — -0.411 (-1.138,
0.390) $\lambda$ 0.261 0.314 0.415 0.482 0.467 0.479 0.476 0.484 0.480 0.476
(0.079, 0.868) (0.095, 1.030) (0.140, 1.312) (0.164, 1.505) (0.160, 1.508)
(0.166, 1.524) (0.159, 1.499) (0.163, 1.565) (0.159, 1.536) (0.153, 1.512)
$\nu$ 1.278 1.277 1.272 1.272 1.273 1.273 1.273 1.272 1.273 1.273 (1.195,
1.363) (1.194, 1.361) (1.187, 1.356) (1.191, 1.357) (1.192, 1.358) (1.191,
1.357) (1.191, 1.357) (1.191, 1.356) (1.190, 1.358) (1.192, 1.356) $\omega$
0.268 0.269 0.273 0.274 0.273 0.273 0.273 0.274 0.273 0.273 (0.203, 0.334)
(0.203, 0.336) (0.209, 0.339) (0.209, 0.339) (0.209, 0.340) (0.210, 0.340)
(0.209, 0.339) (0.209, 0.339) (0.208, 0.341) (0.208, 0.340) $\eta$ -0.113
-0.129 -0.154 -0.168 -0.165 -0.167 -0.166 -0.168 -0.167 -0.166 (-0.219,
-0.009) (-0.234, -0.024) (-0.257, -0.059) (-0.270, -0.073) (-0.269, -0.071)
(-0.270, -0.073) (-0.268, -0.070) (-0.273, -0.073) (-0.271, -0.070) (-0.270,
-0.068)

Table 2: Posterior summaries of the joint model parameters fitted to the NI
renal dataset. Values in the upper rows for each parameter are the 50%
percentiles, whereas the brackets in the lower rows are the 2.5% and 97.5%
percentiles, respectively.

Table 3: Posterior summaries of the joint model parameters fitted to the PBC
dataset. Values in the upper rows for each parameter are the 50% percentiles,
whereas the brackets in the lower rows are the 2.5% and 97.5% percentiles,
respectively.

| nor-nor | mod1 | nor-t-mod2 | nor-t-mod3 | nor-t-tv-1knot | nor-t-tv-2knots | nor-t-tv-3knots | nor-t-tv-4knots | nor-t-tv-5knots
---|---|---|---|---|---|---|---|---|---
$\alpha_{1}$ | 7.851 | 7.468 | 7.479 | 7.515 | 7.521 | 7.517 | 7.521 | 7.526 | 7.529
(7.259, 7.907) | (7.159, 7.779) | (7.161, 7.792) | (7.205, 7.840) | (7.191, 7.848) | (7.194, 7.832) | (7.206, 7.856) | (7.210, 7.863) | (7.197, 7.842)
$\alpha_{2}$ | -0.008 | -0.007 | -0.006 | -0.007 | -0.007 | -0.007 | -0.007 | -0.007 | -0.007
(-0.014, -0.001) | (-0.013, -0.001) | (-0.012, 0.001) | (-0.013, -0.001) | (-0.014, -0.0003) | (-0.013, -0.001) | (-0.014, -0.001) | (-0.014, -0.001) | (-0.013, -0.001)
$\alpha_{3}$ | -0.046 | -0.041 | -0.040 | -0.035 | -0.035 | -0.035 | -0.035 | -0.035 | -0.035
(-0.058, -0.034) | (-0.051, -0.030) | (-0.052, -0.028) | (-0.046, -0.024) | (-0.046, -0.024) | (-0.046, -0.023) | (-0.046, -0.024) | (-0.046, -0.023) | (-0.046, -0.023)
$\alpha_{4}$ | -0.066 | -0.020 | -0.059 | -0.072 | -0.073 | -0.071 | -0.072 | -0.070 | -0.073
(-0.185, 0.054) | (-0.133, 0.095) | (-0.184, 0.062) | (-0.199, 0.051) | (-0.200, 0.049) | (-0.193, 0.054) | (-0.198, 0.050) | (-0.192, 0.051) | (-0.192, 0.046)
$\Sigma_{11}$ | 0.339 | 0.230 | 0.320 | 0.349 | 0.357 | 0.359 | 0.360 | 0.359 | 0.358
(0.286, 0.406) | (0.186, 0.284) | (0.268, 0.384) | (0.295, 0.414) | (0.301, 0.425) | (0.304, 0.426) | (0.303, 0.428) | (0.303, 0.428) | (0.302, 0.427)
$\Sigma_{12}$ | -0.018 | -0.009 | -0.012 | -0.015 | -0.016 | -0.016 | -0.016 | -0.016 | -0.016
(-0.026, -0.011) | (-0.015, -0.004) | (-0.019, -0.005) | (-0.022, -0.009) | (-0.023, -0.009) | (-0.023, -0.009) | (-0.024, -0.010) | (-0.023, -0.010) | (-0.024, -0.010)
$\Sigma_{22}$ | 0.004 | 0.003 | 0.004 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005
(0.003, 0.005) | (0.002, 0.004) | (0.003, 0.006) | (0.004, 0.006) | (0.004, 0.006) | (0.004, 0.006) | (0.004, 0.006) | (0.004, 0.006) | (0.004, 0.006)
$\sigma^{2}$ | 0.107 | 0.047 | 0.040 | 0.021 | 0.023 | 0.024 | 0.024 | 0.024 | 0.024
(0.099, 0.115) | (0.040, 0.054) | (0.034, 0.047) | (0.018, 0.025) | (0.020, 0.027) | (0.021, 0.028) | (0.021, 0.028) | (0.021, 0.029) | (0.021, 0.029)
$\gamma$ | — | 3.516 | — | — | — | — | — | — | —
(2.810, 4.394)
$\delta$ | — | — | 2.759 | 2.091 | — | — | — | — | —
(2.207, 3.451) | (1.824, 2.432)
$\delta_{0}$ | — | — | — | — | 1.449 | 1.184 | 1.130 | 1.122 | 1.130
(1.235, 1.704) | (0.994, 1.410) | (0.940, 1.361) | (0.931, 1.360) | (0.935, 1.367)
$\beta_{1}$ | — | — | — | — | 2.001 | -0.527 | 1.358 | 1.308 | 1.432
(1.432, 2.625) | (-1.408, 0.339) | (0.782, 2.009) | (0.591, 2.118) | (0.664, 2.035)
$\beta_{2}$ | — | — | — | — | -0.543 | 2.946 | 0.202 | 1.376 | 1.175
(-1.459, 0.459) | (2.076, 3.969) | (-0.673, 1.031) | (0.699, 2.128) | (0.391, 2.032)
$\beta_{3}$ | — | — | — | — | — | 1.380 | 2.566 | 0.242 | 1.329
(0.079, 3.056) | (1.720, 3.549) | (-0.816, 1.241) | (0.694, 2.024)
$\beta_{4}$ | — | — | — | — | — | — | 1.210 | 2.414 | 0.314
(-0.192, 3.202) | (1.516, 3.694) | (-0.703, 1.352)
$\beta_{5}$ | — | — | — | — | — | — | — | 1.763 | 2.263
(0.107, 4.256) | (1.321, 3.556)
$\beta_{6}$ | — | — | — | — | — | — | — | — | 1.909
(0.161, 4.574)
$\log(\lambda)$ | -11.276 | -11.427 | -12.129 | -11.952 | -11.901 | -11.911 | -11.861 | -11.854 | -11.872
(-14.850, -7.844) | (-14.814, -8.168) | (-15.584, -8.869) | (-15.215, -8.761) | (-15.300, -8.641) | (-15.256, -8.773) | (-15.178, -8.694) | (-15.200, -8.618) | (-15.242, -8.679)
$\log(\nu)$ | 0.223 | 0.226 | 0.234 | 0.222 | 0.221 | 0.220 | 0.219 | 0.220 | 0.221
(0.064, 0.373) | (0.068, 0.372) | (0.076, 0.380) | (0.069, 0.366) | (0.065, 0.367) | (0.062, 0.367) | (0.065, 0.364) | (0.062, 0.362) | (0.066, 0.363)
$\omega_{1}$ | 0.051 | 0.052 | 0.053 | 0.053 | 0.053 | 0.053 | 0.053 | 0.053 | 0.053
(0.033, 0.070) | (0.035, 0.070) | (0.036, 0.072) | (0.035, 0.071) | (0.035, 0.072) | (0.036, 0.072) | (0.035, 0.071) | (0.036, 0.071) | (0.035, 0.072)
$\omega_{2}$ | -0.157 | -0.121 | -0.144 | -0.137 | -0.137 | -0.136 | -0.136 | -0.133 | -0.139
(-0.519, 0.206) | (-0.460, 0.229) | (-0.508, 0.212) | (-0.488, 0.217) | (-0.492, 0.220) | (-0.498, 0.228) | (-0.494, 0.215) | (-0.483, 0.217) | (-0.490, 0.215)
$\eta$ | 0.799 | 0.805 | 0.895 | 0.873 | 0.872 | 0.869 | 0.865 | 0.862 | 0.865
(0.391, 1.216) | (0.423, 1.193) | (0.511, 1.293) | (0.502, 1.244) | (0.491, 1.253) | (0.497, 1.255) | (0.487, 1.245) | (0.482, 1.250) | (0.488, 1.248)