Handling Missingness Value on Jointly Measured Time-Course and Time-to-event
Data

Gajendra K. Vishwakarma a, Atanu Bhattacharjeeb, Souvik
Banerjeea111Corresponding Author

a Department of Mathematics & Computing,

Indian Institute of Technology Dhanbad, Dhanbad-826004, India,

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
b Section of Biostatistics, Centre for Cancer Epidemiology,

Tata Memorial Center, Navi Mumbai, India

b Homi Bhabha National Institute, Mumbai, India

E-mail<EMAIL_ADDRESS>
###### Abstract

Joint modeling technique is a recent advancement in effectively analyzing the
longitudinal history of patients with the occurrence of an event of interest
attached to it. This procedure is successfully implemented in biomarker
studies to examine parents with the occurrence of tumor. One of the typical
problem that influences the necessary inference is the presence of missing
values in the longitudinal responses as wel l as in covariates. The occurrence
of missingness is very common due to the dropout of patients from the study.
This article presents an effective and detailed way to handle the missing
values in the covariates and response variable. This study discusses the
effect of different multiple imputation techniques on the inferences of joint
modeling implemented on imputed datasets. A simulation study is carried out to
replicate the complex data structures and conveniently perform our analysis to
show its efficacy in terms of parameter estimation. This analysis is further
illustrated with the longitudinal and survival outcomes of biomarkers’ study
by assessing proper codes in R programming language.

Keywords: Joint Modeling, Longitudinal response, Survival outcome, Missing
Data, Multiple Imputation

Mathematics Subject Classification: 62P10

##  1 Introduction

Longitudinal measurements are receiving much attention in clinical research
where patients are followed up to an event of interest. In these studies, the
progression of disease can be closely monitored for a period of time and
outcomes can be predicted more effcetively by assessing temporal changes in
responses. The analysis of follow-up studies with time-to-event observations
is called Joint modeling of longitudinal and survival data (Tsiatis [1995],
Faucett & Thomas [1996])

This article is motivated by the gene expression data from Gene Expression
Omnibus (GEO) database where longitudinal responses were collected for a
number of liquid markers. The patients were given different treatments and the
responses were collected at different time-points. It is examined whether any
progression of tumor is developed for the patients during or after the
treatment period. The problem with the dataset is there are several missing
values in the longitudinal responses as well as covariates which affect the
required inferences. The covariate missingness is an important aspect in this
study as it differentiate the simple response missingness to a more complex
missingness structure. So, it is necessary to develop a proper statistical
methodology to impute those missing values to gather additional information
from the longitudinal and survival patterns.

Gene expression data, collected from micro-array experiments, are useful to
measure temporal changes in gene expression. These data play a crucial role in
understanding the complex mechanism of genes and its effect on the disease
progression. The tumor-specific gene expression profiling has an important
contribution to biomarker detection (Bhattacharjee et al. [2018]). The
potentiality of early detection of disease or providing valuable information
on disease progression presents the protein molecular biomarkers as a popular
choice. Circulating biomarkers, which are obtained in biological liquid
samples, have several advantages over tissue-based profiling due to their non-
invasive nature and efficiency in longitudinal tracking of patients. However,
the cost of collection and measurement of biomarkers is rising and hence the
loss of significant amount of longitudinal responses is unaffordable and it is
essential to retrieve the information lost based on the available data.

During the study related to such a vulnerable disease, it often occurs that
patients leave the study without any prior information (Hui et al. [2013],
Bell et al. [2013]). As a result, this type of dataset contains missing
observations due to drop-out of patients. Missingness can be divided into
three different forms based on dependency pattern on the dataset, namely,
Missing Completely at Random (MCAR), Missing at Random (MAR) and Missing not
at Random (MNAR). In Bayesian and Likelihood methodology, MCAR and MAR
mechanism are well addressed by introducing ignorable missingness and
analyzing the observed data only (Robins et al. [1995], Hogan et al. [2004],
Tsiatis [2007], Bhattacharjee & Vishwakarma [2019]). The MNAR mechanism comes
under the non-ignorable category which requires more specifications of the
joint distribution of the data and missing structure (Ibrahim & Molenberghs
[2009]).

There are methods that simply ignore the covariates for which some
observations are missing and analyze the remaining data as complete case
analysis (Little & Rubin [2014]). But in this process, the sample size is
reduced and efficiency is lost. Also when the dataset is not large enough but
the proportion of missing values is high, the estimates become biased (Sun et
al. [2006], Rubin [1976]). Some attempts have been made to address issues on
covariates in joint modeling. Lu analyzed a joint modeling strategy for skew
longitudinal survival data with missingness and mismeasured covariates using a
logit model (Lu [2017]). A pseudo-likelihood approach was also implemented for
longitudinal binary data with non-ignorable missing response and covariates
(Parzen et al. [2006]). Here a separate model has been considered for the
missingness in covariates. Also, effort has been made in computing joint model
with covariate subject to a limit of detection Sattar & Sinha [2017]. Work has
been performed for missing and left-censored time-varying covariates in joint
modeling using some modified prior distributional assumption (Chen et al.
[2014]). However, there is a clear scarcity of proper modeling technique and
imputation strategy for missing covariates in joint modeling context and it
needs a detailed analytic and computational effort to show a pathway for the
solution of this problem.

In joint models, a classical linear mixed effect model is popularly used to
specify longitudinal measurements and Cox proportional hazard model for the
survival part. The mixed effect model for longitudinal measurements is
associated inside the survival model with an association parameter, (for
example, see (Wulfson & Tsiatis [1997], Henderson [2000], Rizopoulos [2012])).
These characteristics enable joint modeling to capture the influence of
covariates trajectory on the survival duration of patients. Bayesian and
maximum likelihood techniques are efficiently utilized for parameter
estimation (Armero et al. [2016], Song et al. [2016]). The joint modeling work
was stimulated in problems depending on AIDS (Tsiatis [1995]) and currently it
is becoming increasingly popular in other areas also (Ibrahim et al. [2010]).

There are several R software packages which serve as good tools to fit joint
models for a continuous longitudinal data and a time-to-event process under
maximum likelihood method (Rizopoulos [2010]). Bayesian approach has also been
employed in several packages to estimate the parameters and provide more
accurate results in many instances (Rizopoulos [2016]). Some recent packages
use Monte Carlo Expectation-Maximization algorithm for estimation of
parameters and consider multivariate Gaussian process for joint modeling
(Philipson [2018], Hickey et al. [2018]). The mice (Multivariate Imputation of
Chained Equations) package in R provides a good solution for the complex
incomplete data problems (Buuren & Groothuis-Oudshoorn [2010]). It generates
the missing observations using conditional densities by MCMC procedure and
specifies different imputation models for different covariates. Some other
packages use hot-deck imputation, k-nearest neighbor imputation, regression
imputation etc. Other successful efforts considered expectation-maximization
with bootstrapping (EMB) algorithm to impute missing observations considering
the normality assumption and MAR (missing at random) missingness structure
(Honaker [2011]). Iterative robust model-based imputation approach,
implemented through irmi() function in VIM (Kowarik & Templ [2016]) package,
also provides a good alternative to impute missing observations using robust
methods (Templ et al. [2011]). We use R open-source software for the
simulation of posterior samples for the parameters using Markov Chain Monte
Carlo (MCMC) method.

In section 2, we describe the estimation framework of joint models for
longitudinal and time-to-event data and also explain the imputation strategies
for missing data with covariate missingness. In section 3, a simulation study
is described to explain the mentioned methodology. In section 4, we perform
our analysis on a real dataset of gene expression to show its efficacy. In
section 5 and 6, we discuss some results on the proposed technique and also
add future scopes of this methodology.

##  2 Methodology

###  2.1 Joint Modeling

In this section, we have discussed the model-building techniques to analyze
the longitudinal and survival measurements and the joint modeling method to
obtain the efficient treatment strategy.

Longitudinal data is a collection of repeated observations of a particular
subject over a period of time (Weiss [2005]). The temporal ordering of
longitudinal data is important as measurements that are close to each other
tend to similar than the measurements collected at a distant time-point. So,
the analysis also gets difficult with complex variance structure.

One of the widely used tools for analysis of longitudinal measurements is
mixed effect model. The mixed effect model is defined as,

$\begin{split}Y_{i}(t)=M_{i}(t)+\epsilon_{i}(t)=X_{i}(t)\beta+Z_{i}(t)b_{i}+\epsilon_{i}(t),\\\
~{}~{}~{}b_{i}\sim N(0,D),~{}~{}\epsilon_{i}(t)\sim
N(0,\sigma^{2})\end{split}$ (1)

where $Y_{i}(t)$ denotes the longitudinal measurement of $i^{th}$ patient at
time $t$, $M_{i}(t)$ is the true effect of the data or, the true longitudinal
response which we will model in mixed effect model, $\epsilon_{i}(t)$ is the
error component attached to the data. $X_{i}(t)$ & $Z_{i}(t)$ are the fixed
effects and random effects respectively. Also, $\beta$ & $b_{i}$ are the
coefficients of fixed effect and random effects respectively. We assume the
measurement error component $\epsilon_{i}(t)$ follows a Normal distribution
with mean 0 and variance $\sigma^{2}$ and the error is independent of $b_{i}$.
The estimates of the coefficients can be computed using generalized estimating
equations for unknown correlation structures between outcomes.

Dataset consisting of time to the event of interest is referred to as time-to-
event data. When the event of interest is defined as occurrence of death or
any disease to the patient enlisted in the subject, then this kind of dataset
is referred to as Survival Data. Observations are called censored when the
event of interest does not appear for the patient till the time he is
available in the study. Generally, this kind of data is continuous in nature.

Non-parametric methods like Kaplan-Meier estimator is widely used to estimate
and plot the survival probabilities as a function of time. It can also be used
for comparison of multiple groups of subjects. For survival data, often the
interest lies in the association of hazard function and related covariates. A
popularly used parametric model for survival data is Cox proportional hazard
model which is defined as,

$h_{i}(t|K_{i})=h_{0}(t)exp(K_{i}\gamma)$ (2)

where $h_{i}(t|K_{i})$ is the hazard rate for the $i^{th}$ patient who
experiences the event of interest at time t, $h_{0}(t)$ is the baseline hazard
function and $K_{i}$ is baseline covariates with corresponding regression
coefficient $\gamma$. Some popular choices of baseline hazard function are
piecewise constant function, Weibull risk function and unspecified risk
function. In most applications, one of the main objectives is to estimate the
$\gamma$ which explains how the covariates are related to the hazard function.

In joint modeling, we associate both the mixed effect model of longitudinal
measurements and Cox model of survival outcomes in a single model and estimate
the coefficients (Wulfson & Tsiatis [1997], Tsiatis & Davidian [2001]). The
joint model is explained as,

$\begin{split}h_{i}(t|K_{i},M_{i})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
=\lim_{dt\rightarrow 0}Pr(t\leq T_{i}^{*}<t+dt|T_{i}^{*}\geq
t,K_{i},M_{ij})/dt\end{split}$ (3)

$=h_{0}(t)exp(K_{i}\gamma+M_{ij}\alpha)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(4)

Here $\alpha$ measures the effect of longitudinal process on the time-to-event
measurement. It serves as a link between the risk function of cancer infection
at time $t$ with the underlying longitudinal model. The other parameters have
their usual meanings as explained in the longitudinal and survival process.
One of the advantages of this joint model lies in the fact that the parameters
are estimated simultaneously in this approach.

The log-likelihood of the joint model as described by Rizopoulous (Rizopoulos
[2012]) is as follows,

$log~{}p(T_{i},\delta_{i},y_{i};\theta)=log\int
p(T_{i},\delta_{i},y_{i},b_{i};\theta)db_{i}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(5)

$=log\int
p(T_{i},\delta_{i}|b_{i};\theta_{t})[\prod_{j}p\\{y_{i}|b_{i};\theta_{y}\\}]p(b_{i};\theta_{b})db_{i}$
(6)

$\delta_{i}$ is the censoring indicator where 1 implies censored observation
and 0 implies uncensored observation.
$\theta=(\theta_{b},\theta_{t},\theta_{y})$ denotes full parameter vector
where $\theta_{t}$ is parameters for time-to-event data, $\theta_{y}$ denotes
the parameters of longitudinal observations and $\theta_{b}$ denotes the
parameters for random effect component. The parameters of the joint model are
estimated using the maximum likelihood method and the random effects are
predicted using their conditional expectations given all of the data.
Interpretations of the parameters remain the same as what is explained in
random effect model and survival model. Direct maximization of this model is
impossible to work on due to its complex nature and non-parametric components.
The expectation-maximization(EM) algorithm is one of the most used techniques
for optimization purposes of joint likelihood (Wulfson & Tsiatis [1997],
Henderson [2000]).

In bayesian framework, the parameters of the model are estimated using Markov
Chain Monte Carlo(MCMC) method. The posterior distribution of the joint model
is defined as,

$p(\theta,b_{i}|T_{i},\delta_{i},y_{i})\propto
p(T_{i},\delta_{i}|b_{i},\theta)p(y_{i}|b_{i},\theta)p(b_{i}|\theta)p(b_{i}|\theta)p(\theta)$
(7)

where $\theta$ denotes the vector of all the parameters. So, survival
likelihood for patient $i$ is,

$p(T_{i},\delta_{i}|b_{i},\beta,\theta_{s})=h_{i}\\{T_{i}|M_{i}(T_{i}),\theta_{s}\\}^{\delta_{i}}S_{i}\\{T_{i}|M_{i}(T_{i}),\theta_{s}\\}\\\
=[h_{0}(T_{i}|\gamma_{s})exp\\{\gamma^{T}K_{i}+\alpha
M_{i}(T_{i})\\}]^{\delta_{i}}\times
exp\Big{(}-\int_{0}^{T_{i}}h_{0}(s|\gamma_{s})exp\\{\gamma^{T}K_{i}+\alpha
M_{i}(s)\\}ds\Big{)}$ (8)

Here $\theta_{s}=(\gamma_{s},\gamma,\alpha)$ where $\gamma_{s}$ denotes the
parameters associated with the B-splines for the baseline hazard. An inverse
Wishart prior is considered for the variance-covariance matrix ($D$) and a
gamma prior is used for the variance of errors ($\sigma^{2}$).

The Bayesian approach offers the use of Markov Chain Monte Carlo (MCMC)
sampling algorithm to estimate the posterior distribution of the parameters.
The joint models are analyzed using R programming software.

###  2.2 Missing Data Methodology

In real-life situations, obtaining a dataset without any missing observations
especially for patients with diseases is a difficult job. A subject can
refrain itself to be measured at follow-up times which creates a monotone
missing data pattern in longitudinal study. Also, there are some instances
when the subject is absent at one follow-up time and then present again at
subsequent time-points, resulting in intermittent missingness. It is difficult
to evaluate the likelihood functions with non-monotone missing data patterns
as simple factorization rule cannot be applied here. However, in cases where
missing data pattern is considered to be missing at random (MAR), some
imputation models using the software can be applied to estimate the missing
data. Similarly, in survival studies, datasets regarding the death of a
patient often become unavailable generating censored observations. The
covariates or factors impacting upon the failure time are also not available
always throughout the entire study. The difference of measuring frequencies
and the inability to measure all the covariates at a certain time slot will
also result in incomplete data.

A very efficient way to analyze dataset consists of missing observations is
considering the ignorability condition. This condition comes into the picture
when the missingness pattern is MAR and full-data parameters can be decomposed
into two independent parts: one signifies the full data response model and the
other exhibits missing data mechanism. The ignorability condition does not
interpret to drop the missing observations, rather implies that the
missingness mechanism can be left unspecified to draw posterior inference of
the full data model. When the ignorability condition is believed to be not a
suitable option, we can use a more general approach where missingness depends
on the unobserved part of responses given the observed part and covariates.

In MAR, the missingness depends only on the observed values, not on the
unobserved ones (Templ et al. [2011]). In our study, it is found that the
patients tend to skip the clinical tests due to health reasons whenever their
condition is worsening or the size of the tumor is big enough. Also, it can be
found that some of the patients do not feel obliged to test often when values
are not much significant and importance of further testing is not felt. So,
the missingness depends on the observed value of the biomarker i.e., the
missingness is Missing at Random (MAR).

Four different mechanisms that have been used in this article for multiple
imputations are namely, Expectation-Maximization Bootstrap Algorithm,
Predictive Mean Matching, Classification & Regression Tree and Bayesian Linear
Regression Method. The multiple imputation algorithms with covariate
missingness are discussed here. The details on these methods are explained in
the following works (Buuren & Groothuis-Oudshoorn [2010], Honaker [2011],
Rubin [2004], Doove et al. [2014], Vishwakarma et al. [2016]).

#### Classification and Regression Tree (CART)

Cart method is a popular machine learning algorithm which identifies
predictors and cutpoints in the predictors that are used to split the dataset
into more homogeneous subsamples. This process is repeated to obtain a series
of splits that define a binary tree. The imputation algorithm in CART method
is,

(i)Draw a bootstrap sample of $(\dot{y}_{obs},\dot{X}_{obs})$ of size $n_{1}$
from $(y_{obs},X_{obs})$

(ii) Fit $\dot{y}_{obs}$ by $\dot{X}_{obs}$ by a tree model $f(X)$.

(iii) Predict $n_{0}$ terminal nodes from $f(X_{mis})$.

(iv) For $y_{mis}$, determine in which leaf they will end up according to the
tree formed in setp-(iii).

(v) Randomly select one donor $i_{j}$ of the leaf ended up in step (iv) and
replace the missing value $\dot{y}_{j}$ by $y_{j}$, $j=1,2,...,n_{0}$

(vi) Repeat the steps (ii) to (v) several numbers of times.

(vii) The same procedure is performed for each variable in the dataset
replacing the variable containing missing observations with the other
variables for which the data available. Thus at each iteration, we obtain
imputed values for each variable.

#### Expectation-Maximization Bootstrap (EMB) Algorithm

The EMB algorithm is an efficient method for multiple imputation technique.
Here the objective is to break the distribution of observed data and
missingness into other parts where the missingness likelihood can be explained
by observed data. Under MAR conditions, let $D^{obs}$ denotes the observed
data and $M$ the missingness matrix. So, we get,

$p(D^{obs},M|\theta)=p(M|D^{obs})p(D^{obs}|\theta)$ (9)

where $\theta$ is the complete data parameters. In this imputation technique,
it is assumed that the complete data follow multivariate normal distribution
i.e., $D\sim N_{k}(\mu,\Sigma)$. Now, the parameters will be estimated as

$L(\theta|D^{obs})\propto p(D^{obs}|\theta)=\int p(D|\theta)dD^{mis}$ (10)

The data is first divided in $m$ bootstrapped samples. Then for each sample,
the EM algorithm is performed to produce point estimates of $\mu$ and
$\Sigma$. Thus we impute the missing observations using those estimates and
original sample units. Analysis is performed for each of the samples and in
the end, the results are combined by averaging the estimated values.

#### Predictive Mean Matching

In this method, the following algorithm is maintained to impute the missing
observations.

(i) For cases with no missing data, a linear regression of $y$ on $X$ is
fitted and estimates of the coefficients are obtained.

(ii) A random draw is made from the posterior distribution of coefficients.
Typically this can be a random draw from a multivariate normal distribution.

(iii) Using those random draws of the coefficients, the predicted values of
$y$ are generated.

(iv) For each missing $y$, a set of cases with observed $y$ whose predicted
values are close to the observed ones according to any metric defined is
chosen.

(v) From those cases, we randomly choose one and assign its observed value to
the missing observation.

(vi) Repeat the process to generate multiple datasets.

This purpose of this method is to match cases with missing data to similar
cases found in the observed dataset.

#### Bayesian Linear Regression Method

In this method, the missing observations are imputed using Bayesian linear
regression model. The imputation algorithm is,

(i) Calculate $V=(S+diag(S)\kappa)^{-1}$ with small $\kappa$ where
$S=X_{obs}^{\prime}X_{obs}$

(ii) Calculate the regression weights $\hat{\beta}=VX_{obs}^{\prime}y_{obs}$

(iii) Calculate
$\sigma^{2}=(y_{obs}-X_{obs}\hat{\beta})^{\prime}(y_{obs}-X_{obs}\hat{\beta})/g$
where $g$ is a random variable drawn from $\chi^{2}_{(}n-q)$

(iv) Calculate $\dot{\beta}=\hat{\beta}+\dot{\sigma}\dot{z}V^{1/2}$ where
$V^{1/2}$ is computed using Cholesky decomposition and $\dot{z}$ are $q$
random draws from $N(0,1)$

(v) Calculate $\dot{(y)}=X_{mis}\dot{\beta}+\dot{z_{1}}\dot{\sigma}$ where
$\dot{z_{1}}$ is $n_{0}$ independent variates from $N(0,1)$

##  3 Simulation Study

###  3.1 Design

A simulation study was performed to find out the effect of missing
observations on the inference and estimates of parameters. First, we simulated
a full dataset using fixed estimates of the parameters for both longitudinal
and survival model and also selected an association parameter for the joint
model. The submodel for the longitudinal outcome is defined as,

$y_{i}(t)=\beta_{1}+\beta_{2}t+\beta_{3}*(t\\_stage)_{1i}+\beta_{4}*(t\\_stage)_{2i}+\beta_{5}*(t\\_stage)_{3i}+\beta_{6}*(n\\_stage)_{1i}+\\\
b_{i}+\epsilon_{i}(t)$ (11)

where $\beta$’s define the fixed effect parameters and $b$ defines the subject
specific variations which are assumed to be normally distributed
$N(0,\sigma^{2})$. The variance-covariance matrix $D$ ($=\sigma^{2}$, as $b$
is univariate. D is used to keep consistency with the notation) of the random
effects is left unstructured. The $t\\_stage$ and $n\\_stage$ are two
categorical covariates; $t\\_stage$ has 4 levels namely 0,1,2,3 and
$n\\_stage$ has 2 levels namely 0 and 1. The levels are ordinal in nature and
it is assumed that higher the value, greater the effect of the covariate. The
levels of $t\\_stage$ and $n\\_stage$ are drawn using equal probability for
all levels. While choosing the time component, we kept in mind the treatment
plan for patients in the actual scenario. There are 4 patient visit times or
sample collection times considered for simulation of longitudinal outcomes and
the visit times are chosen from the following distributions:
$Poisson(0.7),N(30,5),N(60,5)$ and $N(90,5)$ respectively.

The survival outcomes and times are generated by,

$h_{i}(t)=h_{0}(t)exp(\gamma_{1}*(t\\_stage)_{1i}+\gamma_{2}*(t\\_stage)_{2i}+\gamma_{3}*(t\\_stage)_{3i}+\alpha
m_{i}(t))$ (12)

$m_{i}(t)$ is considered as the true longitudinal response at time $t$. The
starting parameters are drwan randomly from different distributions to
maintain the independence of the regression coefficients and incorporate the
difference in estimation procedure under different methods. The distributions
considered for different parameters are as follows: $\beta_{1}\sim
U(0,10),\beta_{2}\sim U(0.5,1),\beta_{3}\sim U(0.5,1),\beta_{4}\sim
U(0.7,1.2),\beta_{5}\sim U(0.8,1.5),\beta_{6}\sim U(0.5,1.2),\gamma_{1}\sim
U(-0.5,0.5),\gamma_{2}\sim U(-1,1),\gamma_{3}\sim U(-1,1),\sigma\sim
U(0.5,2)$. The baseline hazard function is generated using Weibull
distribution with shape parameter 5 and the association parameter $\alpha$ is
generated using $U(0.5,1.5)$. $D$ or, Var(b) is generated using $U(0.5,2)$.
The censoring time is simulated from an Exponential distribution with mean 2.
To include MAR missingness in the simulation study, the biomarker values are
missing with higher probability for higher the values of $t\\_stage$ and
$n\\_stage$. However, the missingness in covariates is randomly generated. The
covariate missingness is considered to be around 20% of the total
observations.

###  3.2 Analysis

Multiple Imputation (MI) technique is employed to generate the missing data.
Here we use 4 different multiple imputation methods such as Expectation-
Maximization Bootstrap Algorithm, Predictive Mean Matching, Classification &
Regression Tree and Bayesian Linear Regression Method. Imputed datasets are
obtained by using of MICE and Amelia package in R open-source software. For
each of the MI techniques, we obtained 4 imputed datasets and analysis was
carried out on each dataset. Once the entire analysis is carried out, the
results are combined to obtain the estimates of the parameters and they are
compared with the true parameter values considered for simulation.

We genearted 5000 sets of starting values of the paramters and repeated the
simulation study and corresponding analysis. For each generated dataset
containing missing observations, we employed the mentioned missing data
techniques and obtained 200 imputed datasets for each of them to incorporate
the sampling fluctuations among the imputed samples. The joint modeling
analysis is performed on each of the those imputed datasets and results are
combined. Thus, as a whole, the simulation study is replicated 5000 times with
200 imputed samples for each generated data for each multiple imputation
method.

A linear mixed effect model is considered for the longitudinal process with
$t\\_stage$ and $n\\_stage$ are considered as covariates. Proportional hazard
model is considered for the survival outcome with the baseline hazard function
is assumed to follow Weibull distribution. We considered two different
approaches for estimating the parameters of joint models namely maximum
likelihood approach and bayesian approach. Results for both the methods are
compared and the estimated parameters are checked with the original values. In
ML approach of joint modeling, we considered unspecified model for baseline
hazard function whereas for bayesian approach, the baseline risk function is
approximated using splines. Gauss-Hermite quadrature integration was used for
the estimation of parameters in ML approach and Markov Chain Monte Carlo was
used for parameter estimate in bayesian procedure.

The effect of missing observations in response and covariates has been
assessed using two different joint modeling approach. The difference in
parameter values in the simulation study has been examined and on the basis of
the estimates, the multiple imputation technique are judged. The estimates of
the parameters are obtained for each generated dataset and the parameter
estimates (for ML methods) or posterior means (for Bayesian method) and
standard errors (for ML method) or standard deviations (for Bayesian method)
are averaged. As the starting parameters are randomly generated from different
distributions, we compared the performance of different imputation techniques
through two performance measures namely Coverage Probability and Root Mean
Squared Error (RMSE). The performance measures of the parameters for both
response and covariate missingness is displayed in Tables 1 & 2 and the same
for only response missingness are presented in Tables 3 & 4.

Among different multiple imputation techniques, the data imputed by Predictive
Mean Matching (PMM) method results in highest coverage and lowest RMSE among
the four methods followed by EMB algorithm. The Bayesian Normal Regression
(NORM) and Classification and Regression Tree (CART) method performs close to
each other while CART methods performs poorest among those four methods. The
coverages are higher for joint modeling analysis using maximum likelihood
method than bayesian method. Also, the models show better performance when the
missingness is occurred in response only. The RMSEs are much lower for the
association parameter $\alpha$ resulting in better estimate. Among different
imputation methods, most of the estimates obtained from imputed data using EMB
algorithm are found to be close to the full data analysis. For response
missingness only, all the four methods perform better and they are close to
each other while PMM method also performs slightly better than the others.

##  4 Analysis on Gene Expression Study

To perform our analysis, we have used a published dataset on Gene Expression
study collected from Gene Expression Omnibus (GEO) database
(https://www.ncbi.nlm.nih.gov.in/geo/) under the accession number GSE65622. In
this study, all the patients were given neoasjuvant chemotherapy (NACT)
followed by chemoradiotherapy (CRT) and treatment was evaluated for four weeks
after CRT completion. After 2-4 weeks, radical pelvic surgery was considered.
The study consists dataset on 80 patients and the observations were obtained
in the following manner. Serum samples were collected at baseline (Sapling
point=1), following NACT post-NACT(Sampling point=2), at CRT completion
(Sampling point=3), post-CRT (Sampling point=4) and at treatment evaluation
(Sampling point=5). With a high-density antibody array (AHH-BLG-1;RayBiotech
Inc.), the samples were analyzed for the presence of 507 proteins. Further
explanation about the study setting and detailed description of the data can
be cited from Meltzer et al. (2016). The objective of the study is to classify
the proteins according to their effect on cancer progression and investigate
the association between factors regarding tumor severity and duration of
disease reoccurrence.

In the biomarker study, measurements on several biomarkers are gathered from
the patients. It is our interest to find how the missing observations in
$Rantes$, a well-established biomarker related to pulmonary tuberculosis and
HIV disease, affects the analysis of joint modeling. The $Rantes$ biomarker is
also known as $CCL5$ and regulated upon activation, normal T-cell expressed
and secreted. It plays a crucial role in the stimulation of T-cell
proliferation and activates the release of proteins in anti-mycobacterial
immunity Bacon et al. [1995]. The distribution of Rantes is positively skewed,
so a log transformation is applied to the response.The transformed data is
ranging from 3.241 to 5.465 with mean 4.348 and standard deviation 0.326. 65
observations were missing out of 285 records i.e., around 22.81% observations
were missing in the response which consists of 38 patients. Further analysis
shows that the missingness has occurred for 15 out of 28 patients who
experienced the event whereas the responses of 23 out of 51 patients who did
not experience the event are missing. The information regarding the trg_score
are also missing for 3 patients which encourages our study of covariate
missingness in biomarker data. For most of the patients with missing responses
have high t-stage and n-stage values. Hence, we assume that the missingness
structure is dependent on the observed data (MAR missingness). In follow-up
studies of liquid biopsies, patients generally provide responses at early
stages which lead to low number of missing data at the early stages. However,
when they do not find significant improvements in the biomarker values, they
tend to leave the study which leads to missingness. The trajectories of
longitudinal response in shown in the profile plot displayed in Figure 1. Four
different multiple imputation algorithms, as described before, are used for
imputing datasets. 5 datasets are generated for each of the imputation methods
with necessary iteration and the joint modeling is applied on the imputed
dataset. The density comparisons of the imputed data with original data are
shown in Figure 2.

Joint modeling approach was followed to which we explicitly used a linear
mixed effect model for the longitudinal values of the logarithm of biomarker
values. In this study, we considered the longitudinal responses till the
occurrence of event or censoring time, whichever is earlier and rest of the
responses are ignored. In particular, due to randomization technique used in
the study, we considered the fixed effects part of the longitudinal sub-model
as the effect of sampling time (expressed in no. of days and scaled in [0,1]
by dividing with 365), t-stage, n-stage and trg-score. The scores are
independent and hence no interaction term among the covariates is considered.
In the random effects part, we consider a patient-specific random intercept
model.

The longitudinal submodel considered for the gene expression data is,

$log(y_{i}(t))=\beta_{1}+\beta_{2}*(days/365)_{i}+\beta_{3}*t\\_stage_{i}+\beta_{4}*n\\_stage_{i}+\beta_{5}*trg\\_score_{i}+b_{i}+\epsilon_{i}(t)$
(13)

where t_stage, n_stage and trg_score are the ordinal covariates where higher
value denotes more severe effect.

In the survival model, we used Cox Proportional Hazard Model using t-stage and
n-stage as covariates. The survival submodel is defined as,

$h_{i}(t)=h_{0}(t)exp(\gamma_{1}*t\\_stage_{3}+\gamma_{2}*t\\_stage_{4}++\gamma_{3}*n\\_stage_{1}+\gamma_{4}*n\\_stage_{2}+\alpha
m_{i}(t))$ (14)

The baseline hazard maximum likelihood approach is left unstructured. In
bayesian approach, it is modeled with cubic B-splines with 5 knots are being
on the percentiles of the observed event times.

Both bayesian and ML (maximum likelihood) framework has been used for joint
modeling. In ML framework, the baseline risk function is left unspecified
whereas in bayesian framework, the baseline risk function is approximated
using splines. Also in bayesian framework, a simple association structure has
been used where the longitudinal model is associated with the survival
covariates with an association parameter. We have considered only the
observations which have occurred before the event. The estimates of the
parameters are obtained through joint modeling for each generated dataset and
the parameter estimates and standard errors (for ML estimate) or standard
deviations (for bayesian method) are averaged. Comparison for the bayesian
models and ML models was done using Deviance information criterion (DIC) &
Akaike Information Criterion (AIC) values respectively and they are shown in
Tables 5 and 6.

From the tables, we can find that all the imputation models are pretty close
to each other and estimates of some of the parameters are different in
complete case analysis than imputed data analysis. In complete case analysis,
we choose only those patients where all the longitudinal outcomes till the
time of event are available. By comparing the AIC values, it can be said that
the PMM, CART and EMB algorithm imputation methods are very close to each
other whereas NORM method has higher AIC showing a less preferable method. The
association parameter between longitudinal biomarker and survival event
obtained through ML estimate is nearly 1.65 for all the cases except EMB
algorithm for which the estimate is 2.015. However, for bayesian method, it is
approximately 1.10 for all the methods except EMB algorithm for which the
estimate is 1.916.

##  5 Discussion

Joint modeling of longitudinal and survival data with missing covariates is a
challenging area in oncology research. The presence of missing observations
makes it difficult to analyze the data and often leads to biased inference.
Available techniques based on complete observations often leads to biased
estimate when the proportion of missingness is high enough in the dataset. The
joint models combine the conventional linear mixed effect model for
longitudinal responses with the renowned Cox proportional hazard model for the
survival process. Maximizing the likelihood is not always possible in
straightway due to the non-explicit form of the likelihood function which
leads to the use of numerical techniques and bayesian methodology.

It is more challenging to compute the conditional distribution of the
missingness mechanism present in the joint models. Listwise deletion method to
get the complete data is valid only under MCAR conditions, but not for MAR and
MNAR situations. More specifically, when the probability of missingness is
dependent on the unobserved responses and covariates, it requires special
techniques to handle such scenarios.

Separate modeling mechanism is necessary to postulate the cases to obtain
valid inference for joint models. Three main modeling frameworks for missing
observations are a pattern mixture model, shared parameter model, and
selection model. However, these models do not take into account the
missingness of the covariates and are mainly concerned with the response
variable. In our study, we used multiple imputation techniques to obtain the
missing observations in the covariates also. Among several algorithms, we took
support of four methods for data imputation through EMB algorithm, Correlation
and regression tree, Bayesian normal regression, and Predictive mean matching.
The Bayesian perspective in the form of EMB algorithm for imputation has
produced better estimates among those methods. Prior distributional
assumptions are used to generate statistical inferences from the posterior
distributions of the unknown variables including the missing data. Well-known
packages in R programming software were used to analyze the dataset.

##  6 Conclusion

In this article, we demonstrated an effective way to analyze data in joint
modeling scenario when the response and covariates both have missing
observations. Missing covariates for joint longitudinal and survival data
modeling is a challenging task to solve. Our approach shows an efficient way
to handle such missingness in Biomarker study with different multiple
imputation techniques and their effects on joint modeling. Both the joint
modeling approach is assessed to extensively examine the effect of missing
observations on different estimation procedures. In Biomarker study, the cost
of obtaining dataset is high, so it is necessary to use proper techniques to
handle missing observations.

Our procedure is still limited in handling missingness in the parametric
approach. However, further research work is required to accommodate the non-
gaussian response process. The semi-parametric approach could be an attractive
alternative choice for multiple imputations. It is more flexible than others
with extension for non-Gaussian response structure to fit joint models.

##  7 Supplementary Material

The detailed R programming code used in this article is available in the
following link https://github.com/souvikbanerjee91/missing-data-in-joint-
modeling.git

##  8 Acknowledgment

Authors are deeply indebted to the editor Prof. Narayanaswamy Balakrishnan and
learned referee for their valuable suggestions leading to improving the
quality of contents and presentation of the original manuscript. Authors are
also thankful to Science and Engineering Research Board, Department of Science
& Technology, Government of India, for providing necessary support to carry
out the present research work through project Grant No. EMR/2016/003305.

## References

  * Armero et al. [2016] Armero, C., Forné, C., Rué, M., Forte, A., Perpiñán, H., Gómez, G., & Baré, M. (2016) Bayesian joint ordinal and survival modeling for breast cancer risk assessment. Statistics in medicine, 35(28), 5267-5282.
  * Bacon et al. [1995] Bacon, K. B., Premack, B. A., Gardner, P., & Schall, T. J. (1995) Activation of dual T cell signaling pathways by the chemokine RANTES. Science, 269(5231), 1727-1730.
  * Bell et al. [2013] Bell, M. L., Kenward, M. G., Fairclough, D. L., & Horton, N. J. (2013) Differential dropout and bias in randomised controlled trials: when it matters and when it may not. BMJ (Clinical research ed.), 346, e8668.
  * Bernard & Wittwer [2002] Bernard, P. S., & Wittwer, C. T. (2002) Real-time PCR technology for cancer diagnostics. Clinical Chemistry, 48(8), 1178-1185.
  * Bhattacharjee & Vishwakarma [2019] Bhattacharjee, A. and Vishwakarma, G.K., Time-course data prediction for repeatedly measured gene expression data. International Journal of Biomathematics, 12(4), 1950033.
  * Bhattacharjee et al. [2018] Bhattacharjee, A., Vishwakarma, G.K., and Thomas, A., Bayesian state-space modeling in gene expression data analysis: an application with biomarker prediction. Mathematical Biosciences, 305, 96-101, 2018.
  * Buuren & Groothuis-Oudshoorn [2010] Buuren, S. V., & Groothuis-Oudshoorn, K. (2010) mice: Multivariate Imputation by Chained Equations in R. Journal of statistical software, 45(3), 1-67.
  * Chen et al. [2014] Chen, Q., May, R. C., Ibrahim, J. G., Chu, H., & Cole, S. R. (2014) Joint modeling of longitudinal and survival data with missing and left‐censored time‐varying covariates. Statistics in medicine, 33(26), 4560-4576.
  * Cox [1972] Cox, D. R. (1972) Regression models and life‐tables. Journal of the Royal Statistical Society: Series B (Methodological), 34(2), 187-202.
  * Cox & Oakes [1984] Cox, D. R., & Oakes, D. (1984) Analysis of survival data. Chapman & Hall/CRC, London.
  * Doove et al. [2014] Doove, L. L., Buuren, S. V., & Dusseldorp, E. (2014) Recursive partitioning for missing data imputation in the presence of interaction effects. Computational Statistics & Data Analysis, 72, 92-104.
  * Faucett & Thomas [1996] Faucett, C.L. & Thomas, D.C. (1996) Simultaneously modeling censored survival data and repeatedly measured covariates: a Gibbs sampling approach, Statistics in Medicine, 15(15), 1663-1685.
  * Follmann & Wu [1995] Follmann, D. & Wu, M. C. (1995) An approximate generalized linear model with random effects for informative missing data. Biometrics, 51(1), 151–168.
  * Heckman [1976] Heckman, J. J. (1976) The common structure of statistical models of truncation, sample selection and limited dependent variables and a simple estimator for such models. Annals of Economic and Social Measurement, 5(4), 475–492.
  * Henderson [2000] Henderson, R., Diggle, P., Dobson, A. (2000) Joint modeling of longitudinal measurements and event time data., Biostatistics, 1(4),465– 480.
  * Hickey et al. [2018] Hickey, G. L., Philipson, P., Jorgensen, A., & Kolamunnage-Dona, R. (2018) joineRML: a joint model and software package for time-to-event and multivariate longitudinal outcomes. BMC medical research methodology, 18(1), 50.
  * Hogan et al. [2004] Hogan, J. W., Roy, J. and Korkontzelou, C. (2004) Handling drop‐out in longitudinal studies. Statistics in medicine, 23(9), 1455-1497.
  * Honaker [2011] Honaker, J., King, G., & Blackwell, M. (2011) Amelia II: A program for missing data. Journal of statistical software, 45(7), 1-47.
  * Hui et al. [2013] Hui, D., Glitza, I., Chisholm, G., Yennu, S., & Bruera, E. (2013) Attrition rates, reasons, and predictive factors in supportive care and palliative oncology clinical trials. Cancer, 119(5), 1098-1105.
  * Ibrahim & Molenberghs [2009] Ibrahim, J.G., Molenberghs G. (2009) Missing data methods in longitudinal studies: a review. Test, 18(1), 1-43.
  * Ibrahim et al. [2010] Ibrahim J.G., Chu H & Chen L.M. (2010) Basic concepts and methods for joint models of longitudinal and survival data. Journal of Clinical Oncology, 28(16), 2796–2801.
  * Iqbal & Iqbal [2014] Iqbal, N., & Iqbal, N. (2014) Imatinib: A Breakthrough of Targeted Therapy in Cancer, Chemotherapy Research and Practice, Article ID 357027, 9 pages.
  * Kowarik & Templ [2016] Kowarik, A., & Templ, M. (2016) Imputation with the R Package VIM. Journal of Statistical Software, 74(7), 1-16.
  * Little [1994] Little, R. J. A. (1994) A class of pattern-mixture models for normal incomplete data. Biometrika, 81(3), 471–483.
  * Little & Rubin [2014] Little, R. J., & Rubin, D. B. (2014) Statistical analysis with missing data (Vol. 333). John Wiley & Sons.
  * Lu [2017] Lu, T. (2017) Jointly modeling skew longitudinal survival data with missingness and mismeasured covariates. Journal of Applied Statistics, 44(13), 2354-2367.
  * Malvezzi et al. [2017] Malvezzi, M., Carioli, G., Bertuccio, P., Boffetta, P., Levi, F., La Vecchia, C., & Negri, E. (2017) European cancer mortality predictions for the year 2017, with focus on lung cancer. Annals of Oncology, 28(5), 1117-1123.
  * Marrugo-Ramirez et al. [2018] Marrugo-Ramírez, J., Mir, M., & Samitier, J. (2018) Blood-based cancer biomarkers in liquid biopsy: A promising non-invasive alternative to tissue biopsy. International Journal of Molecular Sciences, 19(10), 2877.
  * Molenberghs [2007] Molenberghs, G. & Kenward, M. G. (2007) Missing Data in Clinical Studies (Vol. 61). John Wiley & Sons.
  * Padma [2015] Padma V. V. (2015) An overview of targeted cancer therapy. BioMedicine, 5(4), 19.
  * Parzen et al. [2006] Parzen, M., Lipsitz, S. R., Fitzmaurice, G. M., Ibrahim, J. G., & Troxel, A. (2006) Pseudo‐likelihood methods for longitudinal binary data with non‐ignorable missing responses and covariates. Statistics in medicine, 25(16), 2784-2796.
  * Perkel [2013] Perkel, J. M. (2013) PCR Past, Present & Future. The Scientist Magazine, 27(12), 63-65.
  * Philipson [2018] Philipson P, Sousa I, Diggle P, Williamson P, Kolamunnage-Dona R, Henderson R, Hickey G (2018) joineR: Joint modelling of repeated measurements and time-to-event data. R package version 1.2.4, https://github.com/graemeleehickey/joineR/.
  * Pray [2008] Pray, L. A. (2008) Gleevec: The breakthrough in cancer treatment. Nature Education, 1(1), 37.
  * Rizopoulos [2010] Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software (Online), 35(9), 1-33.
  * Rizopoulos [2012] Rizopoulos D. (2012) Joint models for longitudinal and time-to-event data: With applications in R. Chapman and Hall/CRC, New York.
  * Rizopoulos [2016] Rizopoulos, D. (2016) The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software, 72(7), 1-46.
  * Robins et al. [1995] Robins, J. M., Rotnitzky, A., & Zhao, L. P. (1995) Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association, 90(429), 106-121.
  * Rubin [1976] Rubin, D. B. (1976) Inference and missing data. Biometrika, 63(3), 581-592.
  * Rubin [2004] Rubin, D. B. (2004). Multiple imputation for nonresponse in surveys (Vol. 81). John Wiley & Sons.
  * Sattar & Sinha [2017] Sattar, A., & Sinha, S. K. (2019). Joint modeling of longitudinal and survival data with a covariate subject to a limit of detection. Statistical methods in medical research, 28(2), 486-502.
  * Shih et al. [2015] Shih, Y. C. T., Smieliauskas, F., Geynisman, D. M., Kelly, R. J., & Smith, T. J. (2015) Trends in the cost and use of targeted cancer therapies for the privately insured nonelderly: 2001 to 2011. Journal of Clinical Oncology, 33(19), 2190.
  * Song et al. [2016] Song, H., Peng, Y., & Tu, D. (2016) Recent Development in the Joint Modeling of Longitudinal Quality of Life Measurements and Survival Data from Cancer Clinical Trials. Advanced Statistical Methods in Data Science, 153-168. Springer, Singapore.
  * Sorokin [2000] Sorokin, P. (2000) Mylotarg™ Approved for Patients With CD33+ Acute Myeloid Leukemia. Clinical Journal of Oncology nursing, 4(6).
  * Sun et al. [2006] Sun, Y., Ma, L., Mathew, J., Wang, W., & Zhang, S. (2006) Mechanical systems hazard estimation using condition monitoring. Mechanical systems and signal processing, 20(5), 1189-1201.
  * Templ et al. [2011] Templ, M., Kowarik, A., & Filzmoser, P. (2011). Iterative stepwise regression imputation using standard and robust methods. Computational Statistics & Data Analysis, 55(10), 2793-2806.
  * Tsiatis [1995] Tsiatis, A. A., Degruttola, V., & Wulfsohn, M. S. (1995) Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS. Journal of the American Statistical Association, 90(429), 27-37.
  * Tsiatis & Davidian [2001] Tsiatis, A. A.,Davidian, M., (2001) A semiparametric estimator for the proportional hazards model with longitudinal covariates measured with error, Biometrika, 88(2), 447-458.
  * Tsiatis [2007] Tsiatis, A. (2007) Semiparametric theory and missing data. Springer Science & Business Media, New York.
  * Vishwakarma et al. [2016] Vishwakarma, G.K., Bhattacharjee, A., Jose, J. and Ramesh, V. (2016) Cancer patients missing pain score information: Application with imputation techniques. Epidemiology, Biostatistics and Public Health, 13(4), e11916: 1-8.
  * Weiss [2005] Weiss, R. E., (2005) Modeling longitudinal data. Springer Science & Business Media, New York.
  * Wu & Caroll [1988] Wu, M. C. & Carroll, R. J. (1988) Estimation and comparison of changes in the presence of informative right censoring by modeling the censoring process. Biometrics, 44(1), 175–188.
  * Wulfson & Tsiatis [1997] Wulfsohn, M.S. & Tsiatis, A.A. (1997) A joint model for survival and longitudinal data measured with error, Biometrics, 53(1), 330-339.

Figure 1: Profile plot of log values of the Rantes Biomarker for different
patients

(a)

(b)

(c)

(d)

Figure 2: Density comparisons of imputed data for 4 different imputation
techniques applied on Biomarker data (a) EMB Algorithm, (b) CART method, (c)
NORM method, (d) PMM Method. 5 imputed datasets per method were generated and
plotted with the observed biomarker values for comparison.

Table 1: Coverage Probability and Root Mean Squared Error of Bayesian Estimates of Coefficients in Joint Models on Simulated Biomarker Data with both response and covariate missingness (The RMSEs are given in the brackets) Variable |  | Full
---
Dataset
| EMB
---
Algorithm
| PMM
---
Method
| NORM
---
Method
| CART
---
Method
$\beta_{1}$ | 0.97 (0.13) | 0.81 (0.22) | 0.91 (0.18) | 0.76 (0.24) | 0.72 (0.26)
$\beta_{2}$ | 0.98 (0.21) | 0.74 (0.29) | 0.82 (0.23) | 0.71 (0.28) | 0.67 (0.32)
$\beta_{3}$ | 0.95 (0.13) | 0.74 (0.23) | 0.86 (0.22) | 0.70 (0.28) | 0.64 (0.39)
$\beta_{4}$ | 0.96 (0.17) | 0.73 (0.25) | 0.84 (0.18) | 0.68 (0.30) | 0.61 (0.28)
$\beta_{5}$ | 0.97 (0.22) | 0.68 (0.23) | 0.83 (0.26) | 0.72 (0.31) | 0.66 (0.34)
$\beta_{6}$ | 0.98 (0.11) | 0.78 (0.19) | 0.81 (0.13) | 0.65 (0.22) | 0.62 (0.24)
$\gamma_{1}$ | 0.97 (0.22) | 0.67 (0.31) | 0.78 (0.29) | 0.67 (0.32) | 0.60 (0.33)
$\gamma_{2}$ | 0.98 (0.24) | 0.67 (0.35) | 0.77 (0.29) | 0.72 (0.37) | 0.62 (0.35)
$\gamma_{3}$ | 0.96 (0.23) | 0.72 (0.35) | 0.83 (0.28) | 0.73 (0.34) | 0.71 (0.37)
$\alpha$ | 0.98 (0.15) | 0.71 (0.24) | 0.84 (0.22) | 0.66 (0.28) | 0.68 (0.31)
D | 0.97 (0.06) | 0.73 (0.20) | 0.75 (0.15) | 0.69 (0.26) | 0.62 (0.28)
$\sigma$ | 0.96 (0.11) | 0.68 (0.19) | 0.72 (0.17) | 0.67 (0.23) | 0.61 (0.26)
Table 2: Coverage Probability and Root Mean Squared Error of Maximum Likelihood Estimates of Coefficients in Joint Models on Simulated Biomarker Data with both response and covariate missingness (The RMSEs are given in the brackets) Variable |  | Full
---
Dataset
| EMB
---
Algorithm
| PMM
---
Method
| NORM
---
Method
| CART
---
Method
$\beta_{1}$ | 0.98 (0.12) | 0.85 (0.20) | 0.90 (0.17) | 0.83 (0.23) | 0.80 (0.25)
$\beta_{2}$ | 0.96 (0.25) | 0.82 (0.30) | 0.87 (0.28) | 0.80 (0.31) | 0.78 (0.30)
$\beta_{3}$ | 0.96 (0.15) | 0.82 (0.22) | 0.80 (0.20) | 0.78 (0.26) | 0.76 (0.29)
$\beta_{4}$ | 0.95 (0.21) | 0.82 (0.29) | 0.85 (0.26) | 0.78 (0.31) | 0.76 (0.31)
$\beta_{5}$ | 0.97 (0.20) | 0.82 (0.29) | 0.86 (0.23) | 0.80 (0.30) | 0.76 (0.30)
$\beta_{6}$ | 0.99 (0.13) | 0.80 (0.17) | 0.86 (0.16) | 0.81 (0.21) | 0.80 (0.25)
$\gamma_{1}$ | 0.97 (0.12) | 0.83 (0.17) | 0.89 (0.15) | 0.81 (0.20) | 0.75 (0.24)
$\gamma_{2}$ | 0.96 (0.13) | 0.84 (0.20) | 0.87 (0.18) | 0.81 (0.24) | 0.77 (0.30)
$\gamma_{3}$ | 0.98 (0.18) | 0.83 (0.22) | 0.87 (0.20) | 0.80 (0.26) | 0.79 (0.28)
$\alpha$ | 0.98 (0.14) | 0.82 (0.19) | 0.88 (0.17) | 0.80 (0.24) | 0.76 (0.27)
D | 0.97 (0.13) | 0.79 (0.20) | 0.83 (0.17) | 0.74 (0.24) | 0.73 (0.23)
$\sigma$ | 0.97 (0.06) | 0.81 (0.11) | 0.86 (0.09) | 0.77 (0.12) | 0.72 (0.15)
Table 3: Coverage Probability and Root Mean Squared Error of Bayesian Estimates of Coefficients in Joint Models on Simulated Biomarker Data with only response missingness (The RMSEs are given in the brackets) Variable |  | Full
---
Dataset
| EMB
---
Algorithm
| PMM
---
Method
| NORM
---
Method
| CART
---
Method
$\beta_{1}$ | 0.97 (0.13) | 0.90 (0.20) | 0.92 (0.17) | 0.87 (0.22) | 0.81 (0.26)
$\beta_{2}$ | 0.98 (0.21) | 0.81 (0.26) | 0.87 (0.24) | 0.83 (0.28) | 0.78 (0.29)
$\beta_{3}$ | 0.95 (0.13) | 0.85 (0.23) | 0.91 (0.221) | 0.81 (0.25) | 0.75 (0.31)
$\beta_{4}$ | 0.96 (0.17) | 0.82 (0.21) | 0.87 (0.17) | 0.80 (0.25) | 0.77 (0.25)
$\beta_{5}$ | 0.97 (0.22) | 0.81 (0.23) | 0.84 (0.25) | 0.80 (0.28) | 0.75 (0.31)
$\beta_{6}$ | 0.98 (0.11) | 0.83 (0.13) | 0.86 (0.12) | 0.80 (0.17) | 0.73 (0.20)
$\gamma_{1}$ | 0.97 (0.22) | 0.82 (0.26) | 0.84 (0.27) | 0.80 (0.30) | 0.74 (0.31)
$\gamma_{2}$ | 0.98 (0.24) | 0.78 (0.29) | 0.85 (0.26) | 0.80 (0.30) | 0.76 (0.32)
$\gamma_{3}$ | 0.96 (0.23) | 0.78 (0.28) | 0.84 (0.26) | 0.77 (0.31) | 0.82 (0.35)
$\alpha$ | 0.98 (0.15) | 0.86 (0.18) | 0.90 (0.16) | 0.82 (0.20) | 0.84 (0.25)
D | 0.98 (0.06) | 0.84 (0.15) | 0.87 (0.10) | 0.82 (0.19) | 0.80 (0.22)
$\sigma$ | 0.96 (0.11) | 0.84 (0.14) | 0.88 (0.14) | 0.82 (0.17) | 0.77 (0.19)
Table 4: Coverage Probability and Root Mean Squared Error of Maximum Likelihood Estimates of Coefficients in Joint Models on Simulated Biomarker Data with only response missingness (The RMSEs are given in the brackets) parameters | Full Data |  | EMB
---
Algorithm
| PMM
---
Method
| NORM
---
Method
| CART
---
Method
$\beta_{1}$ | 0.98 (0.12) | 0.90 (0.17) | 0.93 (0.15) | 0.88 (0.23) | 0.83 (0.25)
$\beta_{2}$ | 0.96 (0.25) | 0.89 (0.31) | 0.90 (0.28) | 0.87 (0.30) | 0.86 (0.30)
$\beta_{3}$ | 0.96 (0.15) | 0.87 (0.20) | 0.89 (0.18) | 0.81 (0.23) | 0.80 (0.25)
$\beta_{4}$ | 0.95 (0.21) | 0.85 (0.27) | 0.88 (0.23) | 0.83 (0.29) | 0.81 (0.30)
$\beta_{5}$ | 0.97 (0.20) | 0.88 (0.27) | 0.91 (0.21) | 0.82 (0.28) | 0.80 (0.29)
$\beta_{6}$ | 0.99 (0.13) | 0.86 (0.16) | 0.90 (0.14) | 0.83 (0.20) | 0.81 (0.23)
$\gamma_{1}$ | 0.97 (0.12) | 0.85 (0.17) | 0.91 (0.14) | 0.87 (0.19) | 0.82 (0.22)
$\gamma_{2}$ | 0.96 (0.13) | 0.87 (0.18) | 0.90 (0.17) | 0.84 (0.22) | 0.82 (0.26)
$\gamma_{3}$ | 0.98 (0.18) | 0.86 (0.21) | 0.89 (0.20) | 0.83 (0.25) | 0.81 (0.26)
$\alpha$ | 0.98 (0.14) | 0.90 (0.16) | 0.93 (0.16) | 0.88 (0.19) | 0.84 (0.22)
$D$ | 0.97 (0.13) | 0.86 (0.19) | 0.89 (0.14) | 0.82 (0.23) | 0.80 (0.23)
$\sigma$ | 0.97 (0.06) | 0.88 (0.08) | 0.92 (0.08) | 0.85 (0.11) | 0.84 (0.13)
Table 5: Comparison of Maximum Likelihood estimates of Coefficients in Joint Models on differently imputed Biomarker Data (The standard errors of the estimates are given in the brackets) variables |  | Complete Data
---
(Excluding the
missing
Observations)
| EMB
---
Algorithm
| PMM
---
Method
| CART
---
Method
| NORM
---
Method
(Intercept) | 4.232 (0.131) | 4.283 (0.126) | 4.222 (0.124) | 4.212 (0.134) | 4.258 (0.135)
(days/365) | 0.449 (0.139) | 0.502 (0.088) | 0.499 (0.094) | 0.488 (0.105) | 0.447 (0.108)
t_stage3 | 0.136 (0.115) | -0.334 (0.119) | 0.041 (0.119) | -0.004 (0.125) | 0.011 (0.126)
t_stage4 | -0.066 (0.124) | -0.119 (0.125) | -0.119 (0.125) | -0.142 (0.133) | -0.128 (0.135)
n_stage1 | 0.150 (0.144) | 0.061 (0.116) | 0.161 (0.115) | 0.184 (0.119) | 0.137 (0.123)
n_stage2 | -0.128 (0.101) | -0.140 (0.086) | -0.116 (0.085) | -0.073 (0.089) | -0.148 (0.090)
trg_score2 | 0.067 (0.096) | 0.110 (0.077) | 0.122 (0.078) | 0.112 (0.083) | 0.118 (0.082)
trg_score3 | 0.047 (0.152) | 0.135 (0.096) | 0.135 (0.098) | 0.115 (0.104) | 0.109 (0.107)
trg_score4 | 0.356 (0.163) | 0.149 (0.123) | -0.293 (0.122) | 0.240 (0.127) | 0.243 (0.129)
trg_score5 | 0.488 (0.218) | 0.496 (0.225) | 0.504 (0.240) | 0.564 (0.234) | 0.561 (0.242)
t_stage3 | 0.733 (1.041) | 0.395 (0.974) | 0.256 (0.962) | 0.390 (0.967) | 0.325 (0.972)
t_stage4 | 1.991 (1.041) | 1.439 (0.985) | 1.315 (0.955) | 1.423 (0.964) | 1.364 (0.973)
n_stage1 | 1.991 (0.933) | 1.126 (0.800) | 0.987 (0.807) | 0.989 (0.807) | 0.129 (0.816)
n_stage2 | 0.784 (0.746) | 0.684 (0.605) | 0.641 (0.598) | 0.582 (0.605) | 0.697 (0.608)
| Association
---
Parameter ($\alpha$)
1.334 (0.240) | 2.015 (0.221) | 1.651 (0.214) | 1.646 (0.216) | 1.679 (0.219)
| Variance of
---
random component
(D)
0.193 (0.214) | 0.208 (0.251) | 0.188 (0.247) | 0.565 (0.248) | 0.600 (0.263)
AIC | 188.441 | 415.523 | 419.878 | 413.261 | 466.682
Table 6: Comparison of Bayesian estimates of Coefficients in Joint Models on differently imputed Biomarker Data (The posterior standard deviations of the estimates are given in the brackets) Variable |  | Complete Data
---
(Excluding the
missing
Observations)
| EMB
---
Algorithm
| PMM
---
Method
| CART
---
Method
| NORM
---
Method
(Intercept) | 4.211 (0.428) | 4.273 (0.293) | 4.227 (0.312) | 4.218 (0.287) | 4.242 (0.314)
(days/365) | 0.536 (0.188) | 0.548 (0.154) | 0.562 (0.145) | 0.536 (0.148) | 0.529 (0.165)
t_stage3 | 0.136 (0.410) | 0.006 (0.297) | 0.041 (0.307) | -0.002 (0.285) | 0.019 (0.310)
t_stage4 | -0.033 (0.432) | -0.083 (0.298) | -0.108 (0.316) | -0.115 (0.291) | -0.102 (0.314)
n_stage1 | 0.176 (0.449) | 0.061 (0.283) | 0.159 (0.293) | 0.183 (0.266) | 0.146 (0.285)
n_stage2 | -0.121 (0.313) | -0.143 (0.209) | -0.142 (0.213) | -0.079 (0.198) | -0.144 (0.214)
trg_score2 | 0.071 (0.323) | 0.096 (0.131) | 0.125 (0.155) | 0.102 (0.100) | 0.115 (0.151)
trg_score3 | 0.045 (0.521) | 0.107 (0.163) | 0.149 (0.152) | 0.092 (0.168) | 0.114 (0.164)
trg_score4 | 0.345 (0.550) | 0.171 (0.239) | 0.232 (0.176) | 0.186 (0.168) | 0.216 (0.185)
trg_score5 | 0.433 (0.793) | 0.482 (0.546) | 0.565 (0.626) | 0.545 (0.579) | 0.544 (0.602)
t_stage3 | 1.312 (1.290) | 0.513 (1.091) | 0.508 (1.037) | 0.461 (1.085) | 0.563 (1.103)
t_stage4 | 2.739 (1.295) | 1.553 (1.133) | 1.567 (1.083) | 1.503 (1.066) | 1.631 (1.117)
n_stage1 | 2.482 (1.106) | 1.195 (0.850) | 1.080 (0.841) | 1.094 (0.856) | 1.167 (0.848)
n_stage2 | 1.156 (0.895) | 0.709 (0.643) | 0.669 (0.642) | 0.630 (0.648) | 0.690 (0.630)
| Association
---
Parameter ($\alpha$)
1.842 (0.687) | 1.916 (0.570) | 1.163 (0.655) | 1.078 (0.515) | 1.156 (0.441)
| Variance of
---
random component
(D)
0.728 (0.215) | 0.596 (0.246) | 0.599 (0.248) | 0.565 (0.241) | 0.600 (0.272)
DIC | 341.982 | 665.964 | 667.765 | 646.927 | 730.203