# Estimation and variable selection in a joint model of survival times and longitudinal outcomes with random effects.

Antoine Caillebotte, Estelle Kuhn, Sarah Lemler

► To cite this version:

Antoine Caillebotte, Estelle Kuhn, Sarah Lemler. Estimation and variable selection in a joint model of survival times and longitudinal outcomes with random effects.. 2023. hal-04145010v3

HAL Id: hal-04145010

https://hal.inrae.fr/hal-04145010v3

Preprint submitted on 23 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Estimation and variable selection in high dimension in a causal joint model of survival times and longitudinal outcomes with random effects.

Antoine Caillebotte *email:caillebotte.antoine@inrae.fr*

Université Paris-Saclay, INRAE, UR MaIAGE, UMR GQE-Moulon, France

Estelle Kuhn *email:estelle.kuhn@inrae.fr*

Université Paris-Saclay, INRAE, UR MaIAGE, France,

Sarah Lemler *email:sarah.lemler@centralesupelec.fr*

Université Paris-Saclay, Laboratoire MICS, France,

August 5, 2025

**Abstract:** We consider a joint survival and mixed-effects model to explain the survival time from longitudinal data and high-dimensional covariates in a population. The longitudinal data is modeled using a non linear mixed-effects model to account for the inter-individual variability in the population. The corresponding regression function serves as a link function incorporated into the survival model. In that way, the longitudinal data is related to the survival time. We consider a Cox model that takes into account both high-dimensional covariates and the link function. There are two main objectives: first, identify the relevant covariates that contribute to explaining survival time, and second, estimate all unknown parameters of the joint model. For the first objective, we consider the estimate defined by maximizing the marginal log-likelihood regularized with a $\ell_1$-penalty term. To tackle the optimization problem, we implement an adaptive stochastic gradient to handle the latent variables of the non linear mixed-effects model associated with a proximal operator to manage the non-differentiability of the penalty. We rely on an eBIC model choice criterion to select an optimal value for the regularization parameter. Once the relevant covariates are selected, we re-estimate the parameters in the reduced model by maximizing the likelihood using an adaptive stochastic gradient descent. We provide relevant simulations that showcase the performance of the proposed variable selection and parameter estimation method in the joint model. We investigate the effect of censoring and of the presence of correlation between the individual parameters in the mixed model.

**Key words:** Joint model, nonlinear mixed-effects model, Cox model, high dimension covariates, penalized estimation, adaptive stochastic gradient, proximal operator
# Estimation and variable selection in a joint model of survival times and longitudinal outcomes with random effects.

Antoine Caillebotte *email:caillebotte.antoine@inrae.fr*

Université Paris-Saclay, INRAE, UR MAIAGE, UMR GQE-Moulon, France

Estelle Kuhn *email:estelle.kuhn@inrae.fr*

Université Paris-Saclay, INRAE, UR MAIAGE, France,

Sarah Lemler *email:sarah.lemler@centralesupelec.fr*

Université Paris-Saclay, Laboratoire MICS, France,

May 23, 2024

**Abstract:** We consider a joint survival and mixed-effects model to explain the survival time from longitudinal data and high-dimensional covariates. The longitudinal data is modeled using a nonlinear mixed-effects model, where the regression function serves as a link function incorporated into a Cox model. In that way, the longitudinal data is related to the survival time. Additionally, the Cox model takes into account the inclusion of high-dimensional covariates. The main objectives of this research are two-fold: first, to identify the relevant covariates that contribute to explaining survival time, and second, to estimate all unknown parameters of the joint model. For that purpose, we consider the estimate defined through maximizing a Lasso penalized likelihood. To tackle the optimization problem, we implement a preconditioned stochastic gradient to handle the latent variables of the nonlinear mixed-effects model associated with a proximal operator to manage the non-differentiability of the penalty. We provide relevant simulations that showcase the performance of the proposed variable selection and the parameter estimation method in the joint modeling of the Cox and logistic model.

**Key words:** Joint model, nonlinear mixed-effects model, Cox model, high dimension covariates, penalized estimation, preconditioned stochastic gradient, proximal operator
# 1 Introduction

A very current issue in many fields is the need for a better understanding of the interactions between dependent dynamic phenomena. For example, in medicine, this may involve the dynamics of a patient's tumors in oncology and the effects of anti-cancer treatments administered to the patient. Another example in plant science is the dynamics of plant development in a plot and the spread of an epidemic disease or pests in that plot. The considered phenomena are often complex, both in terms of their modes of interaction and their temporal and spatial dynamics. Moreover, these phenomena are frequently observed in populations of heterogeneous or structured individuals, such as patients or plants.

Mathematical modeling has proven to be a powerful tool for understanding the interactions between multiple dynamic phenomena. It also allows for considering variability present in the observed population of individuals. Joint modeling of several phenomena has demonstrated its effectiveness in several fields, including medicine, pharmacology, and biology (Kerioui et al., 2022). A particular case of joint models concerns the simultaneous modeling of longitudinal data and survival data observed on the same individual based on a causal relation. In this type of joint model, longitudinal data are often modeled by a mixed-effects model (Pinheiro and Bates, 2000; Davidian and Giltinan, 1995), and survival data by a survival model such as the Cox model (Cox, 1972). The latter allows for modeling the instantaneous risk of the survival variable as a function of the covariates. It is also possible to include longitudinal data modeling in a causal way as a covariate in the Cox model via a link function. Several authors have studied this model (Wulfsohn and Tsiatis, 1997, Rizopoulos, 2012, Król et al., 2016).

Inference by maximum likelihood is not direct due to the presence of latent variables in the mixed-effects model and can be performed using numerical tools as for example Expectation Maximization (EM) like algorithms (Wulfsohn and Tsiatis, 1997, Guedj, Thiébaut, and Commenges, 2011, Rizopoulos, 2010, Desmée et al., 2017). The EM-type algorithms, such as the classical Stochastic Approximation Expectation Maximization (SAEM), are the most classical approaches for inferring parameters in the presence of latent variables. They have been developed for estimation in general latent variable models. They are straightforward to implement in the context of a curved exponential family based on sufficient statistics of the model. Moreover, theoretical convergence results have been established in this context. However, when the model does not belong to the exponential family, which is the case for such causal joint models, the methodology is not generic in practice, and the theoretical results fail. Some exponentialization tricks have then been proposed to face this restrictive assumption of the curved exponential family. For instance, one of them consists of considering some unknown parameters as random population variables into a modified model. However, Debavelaere and Allassonnière, 2021 have shown that, in general, the parameter returned by the SAEM algorithm on the modified model is not a maximum likelihood of the initial model, and they have suggested the use of this exponentialization trick with variances of the new random population variables decreasing as the iterations of the algorithm progress. This approach also has limitations in practice due to complex algorithmic settings and tuning. Besides gradient-based methods are another type of approach, often omitted for estimating parameters in latent models Cappé, Moulines, and Rydén, 2005. Recently, Baey et al., 2023 suggested using a preconditioned stochastic gradient algorithm to deal with parameter estimation in the presence of latent variables in very general latent variables models. This approach is particularly interesting when considering models that do not belong to the exponential family, as for causal joint models since theoretical convergence properties are established under general assumptions. Moreover the authors showed that this algorithm performs well for the nonlinear logistic growth mixed-effects model, which can be used to represent some longitudinal data. Note that Bayesian numerical methods have also been
# 1 Introduction

A very current issue in many fields is the need for a better understanding of the interactions between dependent dynamic phenomena. For example, in medicine, this may involve the dynamics of a patient's tumors in oncology and the effects of anti-cancer treatments administered to the patient. Another example in plant science is the dynamics of plant development in a plot and the spread of an epidemic disease or pests in that plot. The considered phenomena are often complex, both in terms of their modes of interaction and their temporal and spatial dynamics. Moreover, these phenomena are frequently observed in populations of heterogeneous or structured individuals, such as patients or plants.

Mathematical modeling has proven to be a powerful tool for understanding the interactions between multiple dynamic phenomena. It also allows for considering variability present in the observed population of individuals. Joint modeling of several phenomena has demonstrated its effectiveness in several fields, including medicine, pharmacology, and biology ([13]). A particular case of joint models concerns the simultaneous modeling of longitudinal data and survival data observed on the same individual. In this type of joint model, longitudinal data are often modeled by a mixed-effects model ([17, 4]), and survival data by a survival model such as the Cox model ([3]). The latter allows for modeling the instantaneous risk of the survival variable as a function of the covariates. It is also possible to include longitudinal data modeling as a covariate in the Cox model via a linking function. Several authors have studied this model ([26], [20], [14]). Due to the presence of latent variables in the mixed-effects model, inference by maximum likelihood can be made via Expectation Maximization (EM) like algorithms ([26], [11], [19], [7]). The EM-type algorithms, such as the classical Stochastic Approximation Expectation Maximization (SAEM), are the most classical approaches for inferring parameters in the presence of latent variables. They have been developed for estimation in general latent variable models. They are straightforward to implement in the context of a curved exponential family based on sufficient statistics of the model. Moreover, theoretical convergence results have been established in this context. However, when the model does not belong to the exponential family, which is the case in our context, the methodology is not generic in practice, and the theoretical results fail.

Some exponentialization tricks have then been proposed to face this restrictive assumption of the curved exponential family. For instance, one of them consists of considering some unknown parameters as random population variables. However, [5] have shown that, in general, the parameter returned by the SAEM on the modified model is not a maximum likelihood of the initial model, and they have suggested the use of this exponentialization trick with variances of the new random population variables that decrease as the iterations of the algorithm progress. This approach also has limitations in practice due to complex algorithmic settings and tuning. The gradient-based methods are another type of approach, often omitted for estimating parameters in latent models. Recently, [2] suggested using a preconditioned stochastic gradient algorithm to deal with parameter estimation in the presence of latent variables. This approach is particularly interesting when considering a model that does not belong to the exponential family, as is the case for our joint model. [2] showed that this algorithm performs well for the nonlinear logistic growth mixed-effects model, which can be used to represent some longitudinal data. Note that Bayesian numerical methods have also been proposed in parallel to estimate the parameters of a joint model. ([18], [21], [13]).

Besides, in many applications, today's technological means allow for collecting high-dimensional explanatory covariates. These may include, for example, genetic markers or omics data. In addition to the wealth of information provided by these covariates, they also generate difficulties in the statistical analysis of models as it is necessary to adapt statistical and numerical approaches to their high dimensionality. One possible approach is to consider a penalized estimator, such as the Lasso ([12], [27]), and adapted numerical methods, such as stochastic proximal gradient ([1], [9]).

In this paper, we consider a joint model that combines, through a link function, a nonlinear mixed-effects model for longitudinal data and a Cox model for the survival times, including covariates of high dimension. Our work aims to select the relevant variables among the high-dimensional covariates in the Cox model part of the joint model based on the whole dataset and then to estimate the model's unknown parameters. For that purpose, we propose an estimate for model parameters, which include a Lasso penalization for the regression parameter of the Cox model. To calculate this estimate in practice, we develop an algorithm combining a preconditioned stochastic gradient to deal with the latent variables in the joint model out of the exponential family and a proximal gradient to handle the non-differentiability of the Lasso penalty used for variable selection in the Cox model part. The proposed algorithm is easy to implement in general joint
proposed in parallel to estimate the parameters of a joint model. (R. Brown and G. Ibrahim, 2003, Rizopoulos, 2016, Keriou et al., 2022).

Besides, in many applications, today's technological means allow for collecting high-dimensional explanatory covariates. These may include, for example, genetic markers or omics data. In addition to the wealth of information provided by these covariates, they also generate difficulties in the statistical analysis of models as it is necessary to adapt statistical and numerical approaches to their high dimensionality. One possible approach to face this high dimensional setting is to consider a parsimonious assumption, for example by assuming that only a small subset of covariates are relevant to explain the variability observed on the survival times. In practice one can consider a regularized estimator, such as the LASSO one (He et al., 2015, Yi, Tang, and Sun, 2022), and adapted numerical methods, such as proximal algorithms (Achab, 2017, Fort, Ollier, and Samson, 2019).

In this paper, we consider a causal joint model that combines through a link function a Cox model for the survival times and a nonlinear mixed-effects model for longitudinal data, including covariates of high dimension. The individual variability in the population is taking into account both in the survival part through the individual covariates and in the longitudinal part through the individual parameter of the mixed model. Our work aims to capitalize on the population approach to efficiently select the relevant covariates among the high-dimensional covariates involved in the Cox model part of the joint model based on the whole dataset and then to estimate in the reduced model the unknown parameters. For that purpose, we propose a procedure combining a variable selection step based on a $\ell_1$-penalized maximum likelihood estimate for the regression parameter of the Cox model with a standard maximum likelihood estimation step in the reduced model. To calculate the penalized estimate in practice, we develop an algorithm combining an adaptive stochastic gradient descent to deal with the presence of individual parameters in the mixed model part which are latent variables and a weighted proximal algorithm to handle the non-differentiability of the $\ell_1$ penalty term used for variable selection in the Cox model part. The proposed algorithm is easy to implement in general joint models, in particular it is not require that the model density belongs to the curved exponential family.

The paper is organized as follows. In Section 2, we present the causal joint model constructed from a Cox model for survival data, with high-dimensional covariates and a link function related to a non linear mixed-effects model for longitudinal data; we detail some examples. In Section 3, we present the proposed inference method based on: i) a variable selection procedure composed of the evaluation of the $\ell_1$-penalized maximum likelihood estimate for several value of the regularization parameter on a grid and of the choice of this regularization parameter using the eBIC model choice criterion; ii) the re-estimation of the parameters in the reduced model once the relevant covariates have been selected. The fourth section details the numerical steps, in particular the adaptive proximal weighted stochastic gradient algorithm used to target the $\ell_1$-penalized estimate on the grid. Finally, we illustrate through a simulation study the performance of the proposed methodology in Section 5. The paper ends with a conclusion.

# 2 Joint causal model for survival times and longitudinal outcomes

Let us begin by introducing some notation. We consider a population of $N$ individuals. For each individual $i$, we denote by $T_i^*$ the survival time corresponding to the duration until the occurrence of an event of interest. We also observe for each individual $i$ longitudinal data, which are repeated observations over times. We denote by $Y_{i,j}$ the $j$th observation of the $i$th individual for $i \in \{1, \dots, N\}$ and $j \in \{1, \dots, J\}$. Note that we consider in our work the same number of repeated observations per individual; however it can be easily generalized to the case where there are different numbers of
models without assuming that model density belongs to the curved exponential family.

The paper is organized as follows. In Section 2, we detail the joint model constructed from a nonlinear mixed-effects model for longitudinal data and a Cox model for survival data, with high-dimensional covariates and a link function. In Section 3, we present the proposed inference method based on a Lasso penalized estimator and numerical procedure based on a stochastic proximal gradient algorithm. Finally, we illustrate the methodology in Section 4 through a simulation study. The paper ends with a conclusion.

# 2 Joint model for survival times and longitudinal outcomes

We consider $N$ individuals, and for each individual $i$, we denote the survival time $T_i^*$, corresponding to the duration until the occurrence of an event of interest, and longitudinal data, which are repeated observations $J$ times denoted by $Y_{i,j}$ with $i \in \{1, \dots, N\}$ and $j \in \{1, \dots, J\}$. Note that our work can easily be generalized to the case where there are different numbers of longitudinal observations for each individual in the population. The following describes the joint model we considered.

## 2.1 Survival Model

The survival time $T_i^*$ of individual $i$ is the time between a fixed initial moment and the occurrence of an event of interest. To characterize the distribution of $T_i^*$, we use the hazard function defined by:

$$h_i(t) := \lim_{\Delta t \to 0} \frac{\mathbb{P}(t \le T_i^* < t + \Delta t | T_i^* \ge t)}{\Delta t}; \forall t \ge 0. \quad (1)$$

The Cox model ([3]) is one of the most classical models in survival analysis. It allows us to relate the hazard function of the survival time $T_i^*$ to the covariates $U_i \in \mathbb{R}^p$, with $p$ being the number of covariates. In our approach, we will consider the high-dimensional setting with many covariates so that $p$ is very large with respect to $N$. The hazard of the Cox model for individual $i$ is written as follows:

$$h(t|U_i) = h_0(t) \exp(\beta^T U_i); \forall t \ge 0, \quad (2)$$

with $\beta \in \mathbb{R}^p$ a regression parameter and $h_0$ the baseline hazard function that characterizes a common behavior in the observed population. In the sequel, we will consider a parametric baseline function denoted by $h_{\theta_b}$ where $\theta_b$ are its parameters. Therefore, the Cox model's unknown parameters are $\beta$ and $\theta_b$. It is a positive random variable. In some settings the survival time is not observed directly due to censorship. Thus let's suppose that $T_i^*$ is censored on the right, and let's denote by $C_i$ the censoring time. We then observe $T_i = \min(T_i^*, C_i)$ and $\delta_i = \mathbb{1}_{\{T_i^* \le C_i\}}$ respectively the censored survival time and the censoring indicator.

In addition to the covariates, we consider explaining some of the survival time variability using the longitudinal data dynamic, which will be modeled using a nonlinear mixed-effects model. Let us present the mixed-effects model before explaining the integration of this new component into the Cox model.

## 2.2 Nonlinear Mixed-Effects Model

The longitudinal data are observed $J$ times for each individual $i \in \{1, \dots, N\}$. Let us denote by $Y_{i,j}$ the $j$-th observation of the $i$-th individual for $j \in \{1, \dots, J\}$ and $i \in \{1, \dots, N\}$. We model this longitudinal observation using a nonlinear function $m$ that depends on individual parameters represented by the latent variable $Z_i$ as follows:

$$\left\{ \begin{array}{l} \mathbf{Y}_{i,j} = m(t_j; \mathbf{Z}_i) + \varepsilon_{i,j}, \forall 1 \le i \le N, 1 \le j \le J \\ \mathbf{Z}_i \sim \mathcal{N}(\mu, \Gamma); \varepsilon_{i,j} \sim \mathcal{N}(0, \sigma^2), \end{array} \right. \quad (3)$$

where, $t_j$ is the $j$-th observation time, and $\varepsilon_{i,j}$ is an additive noise assumed centered Gaussian with unknown variance $\sigma^2$. The latent variable $Z_i$ describes the inter-individual variability of the population. It is assumed that $Z_i$ follows a Gaussian distribution with unknown expectation $\mu$ and variance $\Gamma$. The unknown parameters of the nonlinear mixed-effects model are therefore $\mu, \Gamma$, and $\sigma^2$.

Let us introduce in the following the link function, which will combine the two previous models by modeling the influence of the dynamic of the longitudinal observation of the hazard function.
longitudinal observations for each individual in the population. The following sections describe the considered joint model.

## 2.1 Survival Model

For $i \in \{1, \dots, N\}$ the survival time $\mathbf{T}_i^*$ of individual $i$ represents the time between a fixed initial moment and the occurrence of an event of interest and is modeled through a random variable. To characterize the distribution of $\mathbf{T}_i^*$, we consider the hazard function defined by:

$$h_i(t) := \lim_{\Delta t \to 0} \frac{\mathbb{P}(t \le \mathbf{T}_i^* < t + \Delta t | \mathbf{T}_i^* \ge t)}{\Delta t}; \forall t \ge 0. \quad (1)$$

The Cox model (see Cox, 1972) is one of the most classical models in survival analysis. It allows us to relate the hazard function of the survival time $\mathbf{T}_i^*$ to covariates $X_i \in \mathbb{R}^p$, with $p$ being the number of covariates, characterizing individual $i$. In our approach, we will consider the high-dimensional setting with many covariates so that $p$ is very large with respect to $N$. The hazard of the Cox model for individual $i$ is written as follows:

$$h(t|X_i) = h_0(t) \exp(\beta^T X_i); \forall t \ge 0, \quad (2)$$

with $\beta \in \mathbb{R}^p$ a regression parameter and $h_0$ the baseline hazard function that characterizes a common behavior in the observed population. In the sequel, we consider a parametric baseline function denoted by $h_{\theta_b}$ where $\theta_b$ are its parameters. Therefore, the Cox model's unknown parameters are $\beta$ and $\theta_b$.

In many settings the survival time is not observed directly due to censorship (see Cox, 1972). Thus let's suppose that $\mathbf{T}_i^*$ is censored on the right, and let's denote by $C_i$ the censoring time. We then observe $\mathbf{T}_i = \min(\mathbf{T}_i^*, C_i)$ and $\delta_i = \mathbb{1}_{\mathbf{T}_i^* \le C_i}$ respectively the censored survival time and the censoring indicator.

An intrinsic assumption inherent to the Cox model is that the hazards are proportional, meaning that the ratio between the hazard of two individuals only depends on a function of the covariates of these two individuals and does not depend on the time. This assumption is quite restrictive, in particular many dataset do not satisfy it in practice. For this reason, several extensions have been developed to relax this assumption. It is possible for example to include time-dependent covariates or to consider a random effect as in the frailty models. When additional longitudinal data related to the survival data are available, another way to release this assumption is to explain some of the survival time variability using the longitudinal data dynamic, by modeling them with a nonlinear mixed-effects model and including a causal effect through a link into the Cox model.

Let us present the mixed-effects model before explaining the integration of this new component into the Cox model.

## 2.2 Non linear mixed-effects model

The longitudinal data are observed $J$ times for each individual $i \in \{1, \dots, N\}$. Let us denote by $\mathbf{Y}_{i,j}$ the $j$th observation of the $i$th individual for $j \in \{1, \dots, J\}$ and $i \in \{1, \dots, N\}$. We model this longitudinal observation using a non linear function $m$ that depends on individual parameters represented by the latent variable $Z_i$ as follows:

$$\left\{
\begin{array}{ll}
\mathbf{Y}_{i,j} = m(t_j; Z_i) + \varepsilon_{i,j}, & \forall 1 \le i \le N, 1 \le j \le J \\
Z_i \sim \mathcal{N}(\mu, \Gamma) & \forall 1 \le i \le N \\
\varepsilon_{i,j} \sim \mathcal{N}(0, \sigma^2) & \forall 1 \le i \le N, 1 \le j \le J
\end{array}
\right.
\quad (3)$$
## 2.3 Joint model

We introduce the link function denoted by $f$, which is linearly parameterized with $\alpha$ and depends on the history of the true unobserved longitudinal data $\mathcal{M}(t; \mathbf{Z}_i) = \{m(s; \mathbf{Z}_i)|\forall s, 0 \le s < t\}$. $\mathcal{M}(t; \mathbf{Z}_i)$ describes the past values of the longitudinal dynamic up to time $t$. Several forms for $f$ can be taken, an example of which is given below. We assume that the hazard of the survival time of individual $i$ is related to the longitudinal data dynamic through the link function $f$ as follows, $\forall 1 \le i \le N, \forall t \ge 0$:

$$h(t|\mathcal{M}(t, \mathbf{Z}_i), U_i) = h_{\theta_b}(t) \exp(\beta^T U_i + f(\alpha, \mathcal{M}(t; \mathbf{Z}_i))), \quad (4)$$

The parameter $\alpha$ represents the influence of the longitudinal dynamic on the survival data.

The unknown parameters for the joint model include the parameters of the Cox model and those of the nonlinear mixed-effects model, as well as the link function parameter of the joint model. We note $\theta = (\theta_b, \beta, \alpha, \mu, \Gamma, \sigma^2) \in \Theta$ the vector of unknown parameters with $\Theta \subset \mathbb{R}^d$ being the parameter space.

# 3 Inference Method

In this section, we propose a method for estimating the model parameters presented above. We will define two estimators to deal first with the latent variable context and then with the high dimension.

## 3.1 Expression of the Marginal Likelihood

We consider the maximum likelihood estimator to infer the joint model parameters. In the context of latent variable models, the marginal likelihood, denoted by $\mathcal{L}_{\text{marg}}$, is obtained by integrating the complete likelihood over the latent variables, which are not observed. Let $\mathbf{D}_i = (\mathbf{Y}_i, \mathbf{T}_i, \delta_i)$ be the observed variables for $1 \le i \le N$ and $\mathbf{D} = (\mathbf{D}_i)_{1\le i \le N}$:

$$\begin{aligned} \mathcal{L}_{\text{marg}}(\theta; \mathbf{D}) &= \prod_{i=1}^{n} \int p_{\theta}(\mathbf{D}_i, \mathbf{Z}_i) d\mathbf{Z}_i \\ &= \prod_{i=1}^{n} \int p_{\theta}(\mathbf{D}_i|\mathbf{Z}_i) p_{\theta}(\mathbf{Z}_i) d\mathbf{Z}_i, \end{aligned} \quad (5)$$

where $p_{\theta}(\mathbf{D}, \mathbf{Z}), p_{\theta}(\mathbf{D}|\mathbf{Z}), p_{\theta}(\mathbf{Z})$ are respectively the density of the pair $(\mathbf{D}, \mathbf{Z})$, the density of $\mathbf{D}$ conditionally to $\mathbf{Z}$, and the density of $\mathbf{Z}$.

$$\begin{aligned} \mathcal{L}_{\text{marg}}(\theta; \mathbf{D}) &= \prod_{i=1}^{n} \int \left( \left\{ \prod_{j=1}^{J} p_{\theta}(Y_{i,j} | \mathbf{Z}_i) \right\} \right. \\ & \qquad \left. \times p_{\theta}(\mathbf{T}_i, \delta_i | \mathcal{M}(t, \mathbf{Z}_i), U_i) p_{\theta}(\mathbf{Z}_i) \right) d\mathbf{Z}_i \\ &= \prod_{i=1}^{n} \int \left( \left\{ \prod_{j=1}^{J} \Phi(Y_{i,j}; m(t_j, \mathbf{Z}_i), \sigma^2) \right\} \right. \\ & \qquad \left. \times (h(T_i | \mathcal{M}(T_j, \mathbf{Z}_i), U_i))^{\delta_i} \right. \\ & \qquad \left. \times \exp \left( - \int_{0}^{T_i} h(s | \mathcal{M}(s, \mathbf{Z}_i), U_i) ds \right) \right. \\ & \qquad \left. \times \Phi(\mathbf{Z}_i; \mu, \Gamma) \right) d\mathbf{Z}_i, \end{aligned}$$

where $\Phi(.; m, v)$ is the density of normal distribution centered in $m$ and of variance $v$.
where, $t_j$ is the jth observation time, and $\varepsilon_{i,j}$ is an additive noise assumed centered Gaussian with unknown variance $\sigma^2$. The latent variable $Z_i$ describes the inter-individual variability of the population. It is assumed that $Z_i$ follows a Gaussian distribution with unknown expectation $\mu$ and covariance matrix $\Gamma$. We assume that for each $i$ the variables $Z_i$ and $(\varepsilon_{i,j})_j$ are independent. We assume also that the variables $(Z_i)_i$ are independent, as well as the variables $(\varepsilon_{i,j})_{i,j}$. The unknown parameters of the non linear mixed effects model are therefore $\mu, \Gamma$ and $\sigma^2$.

Let us introduce in the following the link function, which will combine the two previous models by modeling the influence of the dynamic of the longitudinal observation of the hazard function.

## 2.3 Joint causal model

We introduce the link function denoted by $f$, which is linearly parameterized with $\alpha$ and depends on the history of the true unobserved longitudinal dynamic. Let us define the set $\mathcal{M}(t; Z_i) = \{m(s; Z_i)|\forall s, 0 \le s < t\}$ that contains the values of the longitudinal dynamic up to time $t$. Several forms for $f$ can be taken, an example of which is given below. We assume that the hazard of the survival time of individual $i$ is related to the longitudinal dynamic through the link function $f$ as follows, $\forall 1 \le i \le N, \forall t \ge 0$:

$$h(t|\mathcal{M}(t, Z_i), X_i) = h_{\theta_b}(t) \exp(\beta^T X_i + f(\alpha, \mathcal{M}(t; Z_i))). \quad (4)$$

The parameter $\alpha$ represents the influence of the longitudinal dynamic on the survival data. The unknown parameters for the joint model include the parameters of the Cox model and those of the non linear mixed-effects model, as well as the link function parameter of the joint model. We denote by $\theta = (\theta_b, \beta, \alpha, \mu, \Gamma, \sigma^2)$ the vector of unknown parameters taking value in the parameter space denoted by $\Theta \subset \mathbb{R}^d$.

# 3 Variable selection and estimation procedure

In this section, we propose a method for selecting the relevant covariates and estimating the model parameters.

## 3.1 Inference in latent variables models

We consider the maximum likelihood estimator to infer the joint model parameters. Since the random effects of the mixed model part are non observed, it is not possible to evaluate the full complete likelihood of the joint model depending on both observations and latent variables. Therefore as done usually in the context of latent variable models, we consider the marginal likelihood, denoted by $\mathcal{L}_{\text{marg}}$ obtained by integrating the full complete likelihood over the latent variables. Let $\mathcal{D}_i = (\mathbf{Y}_i, \mathbf{T}_i, \delta_i)$ be the observed variables for $1 \le i \le N$ and $\mathcal{D} = (\mathcal{D}_i)_{1 \le i \le N}$. The marginal likelihood writes:

$$\mathcal{L}_{\text{marg}}(\theta; \mathcal{D}) = \prod_{i=1}^{n} \int p_{\theta}(\mathcal{D}_i, Z_i) dZ_i = \prod_{i=1}^{n} \int p_{\theta}(\mathcal{D}_i|Z_i) p_{\theta}(Z_i) dZ_i,$$

where $p_{\theta}(\mathcal{D}, \mathbf{Z}), p_{\theta}(\mathcal{D}|\mathbf{Z}), p_{\theta}(\mathbf{Z})$ are respectively the density of the pair $(\mathcal{D}, \mathbf{Z})$, the density of $\mathcal{D}$ conditionally to $\mathbf{Z}$, and the density of $\mathbf{Z}$. We assume that the variables $(\mathbf{Y}_i)$ and $(\mathbf{T}_i)$ are independent conditionally to the variables $(\mathbf{Z}_i)$. We assume also that the censoring times $(C_i)$ and $(\mathbf{T}_i)$ are independent conditionally to the variables $(\mathbf{Z}_i)$. Thus we get the following expression for the likelihood:
Due to the integral, it is difficult to directly compute the maximum of the marginal likelihood, which does not have an analytical form in this latent variable model. Therefore, we will use numerical methods to solve this maximization problem.

## 3.2 Definition of the Penalized Estimator for Variable Selection

One can usually estimate the model parameters by maximizing the likelihood using the maximum likelihood estimator written as follows:

$$ \hat{\theta}_{\text{MLE}}(\mathcal{D}) = \arg\max_{\theta \in \Theta} (\log \mathcal{L}_{\text{marg}}(\theta; \mathcal{D})), \quad (6) $$

where $\Theta$ denotes the parameter space. However, in our context we must deal with the high dimension of the covariates, so we introduce a penalty and consider a penalized maximum likelihood estimator. We aim to select relevant variables among the covariates of the survival model. We use the Lasso (Least Absolute Shrinkage and Selection Operator) procedure, which was initially developed for linear regression models ([24]) and the Cox model ([25]). This method enables us to handle high-dimensional data and select a subset of explanatory covariates from a large collection. We consider a Lasso penalty, which only depends on the parameter $\beta$:

$$ \text{pen}_{\lambda}(\theta) = \lambda \|\beta\|_{1} = \lambda \sum_{k=1}^{p} |\beta_{k}|, $$

where $\lambda$ is a positive real called regularization parameter. Our goal is then to maximize the logarithm of the marginal likelihood where the penalty is integrated as follows. Let us define the penalized maximum likelihood estimator by:

$$ \hat{\theta}_{\text{Lasso}}(\mathcal{D}) = \arg\max_{\theta \in \Theta} (\log \mathcal{L}_{\text{marg}}(\theta; \mathcal{D}) - \text{pen}_{\lambda}(\theta)), \quad (7) $$

where $\Theta$ denotes the parameter space and where $\lambda$ is a positive parameter. The larger the value of $\lambda$, the more $\beta$ will be constrained to have zero components. Conversely, the smaller the value of $\lambda$, the freer the components of $\beta$ will be. It is customary to determine the value of $\lambda$ using cross-validation ([24]). But here, without any predictive context, we study the BIC criterion to find the best regularization value ([23]).

Usually, when we deal with latent variables, since the marginal likelihood is non-analytic, classical methods used to infer the unknown parameters are Expectation Maximization like algorithms ([16]). The inconvenience of these procedures is that they are well adapted to models of the curved exponential family, which is not the case for the joint model we consider. The nonlinearity of the function $m$ on the individual parameters $\varphi$ prevents it from being in this exponential family. Recently [2] have presented a preconditioned stochastic gradient descent for estimation in a latent variable model adapted to general latent variables models. Moreover, due to the non-differentiability of the considered penalty, we will use a proximal algorithm as presented by [1] and [9]. Thus, we add a proximal gradient in the procedure presented in [2] and implement a preconditioned stochastic proximal gradient algorithm to calculate the estimator.

## 3.3 Implementation of the Inference Procedure

We deal simultaneously with unobserved random effects of the mixed-effects model and the penalty term by implementing a Preconditioned Stochastic Proximal Gradient, SPG-FIM, in the sequel. The latter is a forward-backward splitting algorithm that can compute the penalized estimate (7). The algorithm is divided into three steps; a realization of the latent variables is sampled with a first step called *Simulation*, which uses a Metropolis-Hastings sampler ([10]). The second step is the classical gradient descent on the approximate complete likelihood, the *Forward* step. Following the procedure presented in [2], we have chosen to use a preconditioning of the gradient with an estimate of the Fisher information matrix. The latter is updated during the iterations using the estimate presented by [6]. The last step, called *Backward*, deals with the penalty term. We apply the classical proximal operator ([15, 22], defined below:
$$
\begin{align*}
\mathcal{L}_{\text{marg}}(\theta; \mathbf{D}) &= \prod_{i=1}^{N} \int p_{\theta}(\mathbf{Y}_i|Z_i)p_{\theta}(\mathbf{T}_i, \delta_i|Z_i)p_{\theta}(Z_i)dZ_i \\
&= \prod_{i=1}^{N} \int \left( \left\{ \prod_{j=1}^{J} \Phi(Y_{i,j}; m(t_j, Z_i), \sigma^2) \right\} \times (h(T_i|\mathcal{M}(T_j, Z_i), X_i))^{\delta_i} \right. \\
&\qquad \left. \times \exp \left( -\int_{0}^{T_i} h(s|\mathcal{M}(s, Z_i), X_i)ds \right) \times \Phi(Z_i; \mu, \Gamma) \right) dZ_i,
\end{align*}
$$

where $\Phi(.; m, v)$ is the density of the normal distribution with expectation $m$ and covariance $v$.

Model parameters are typically estimated using the maximum likelihood estimator, defined as:

$$
\hat{\theta}_{\text{MLE}}(\mathbf{D}) = \arg \max_{\theta \in \Theta} \log \mathcal{L}_{\text{marg}}(\theta; \mathbf{D}), \quad (5)
$$

where $\Theta$ denotes the parameter space.

However, in practice, it is difficult to directly compute this estimate, which does often not
have an analytical form in latent variable models. Therefore, one has to use numerical methods.
When dealing with latent variables, traditional approaches for computing this estimate often rely
on Expectation-Maximization (EM) algorithms (see Ng, Krishnan, and McLachlan, 2012). A key
limitation of these methods is that they are primarily suited to models within the curved exponential
family, both from theoretical and practical points of view. Indeed theoretical results are established
only for models of the curved exponential family. Moreover implementation outside the exponential
family is complex and require to store numerous quantities in memory. To overcome this limitation,
a common strategy involves applying the exponentiation trick, which consists of modeling some
parameters as Gaussian random variables to construct an augmented model belonging to the
exponential family. Nonetheless, this approach also has limitations, as discussed by Debavelaere and
Allassonnière, 2021. Another way to address this issue is to implement stochastic gradient descent
algorithms to solve the optimization problem using the Fisher identity (see Cappé, Moulines, and
Rydén, 2005). In particular a preconditioned stochastic gradient descent algorithm for parameter
estimation in latent variable models was recently proposed by Baey et al., 2023. This method is
more generally applicable and theoretical convergence results are stated without requiring that the
model belongs to the curved exponential family.

3.2 Definition of the penalized estimator for variable selection

In the presence of covariates of high dimension, classical estimation strategies are no more efficient,
so we introduce a penalty term to regularize the estimation and perform the variable selection.
We borrow the idea of the LASSO (Least Absolute Shrinkage and Selection Operator) procedure,
developed initially for linear regression models (see Tibshirani, 1996) and extended in several models,
including the Cox model (see Tibshirani, 1997). This method enables to handle high-dimensional
covariates and to select a subset of explanatory covariates from a large collection. Assuming only a
few variables are relevant to be included in the model among the whole collection, we consider a $\ell_1$
penalized maximum likelihood estimator to select them. We emphasize that we consider a penalty
term, which only depends on the parameter $\beta$:

$$
\text{pen}_{\lambda}(\theta) = \lambda \|\beta\|_{1} = \lambda \sum_{k=1}^{p} |\beta_{k}|,
$$
$$ \mathrm{Prox}_{\mathrm{pen}_{\lambda}}(\beta) = \arg \min_{\beta' \in \mathbb{R}^p} \left( \mathrm{pen}_{\lambda}(\beta') + \frac{1}{2} \|\beta - \beta'\|_2^2 \right). \qquad (8) $$

With the Lasso penalty, the proximal operator has an explicit form:

$$ (\mathrm{Prox}_{\mathrm{pen}_{\lambda}}(\beta))_i = \begin{cases} 0 & \text{if } |\beta_i| < \lambda \\ \beta_i - \lambda & \text{if } \beta_i \ge \lambda \\ \beta_i + \lambda & \text{if } \beta_i \le -\lambda \end{cases} \quad ; \quad \forall i \in \{1, \dots, p\}. \qquad (9) $$

The *Backward* step corresponds to the application of the proximal operator on the result of the *Forward* step.

As the penalty only depends on $\beta$, the proximal operator selects the $\beta$ components that seem to be the most explanatory of the data. It computes a sparse solution for $\beta$ but also applies shrinkage on the non-zero components so that the Lasso estimator is biased. Therefore, we detail a method to obtain an unbiased estimator in what follows.

Algorithm 1 provides the steps of the stochastic proximal gradient, where $(\gamma_k)_{k \ge 1}$ is a step size such that $\forall k \in \mathbb{N}, \gamma_k \in [0, 1], \sum_{k=1}^{\infty} \gamma_k = \infty$ and $\sum_{k=1}^{\infty} \gamma_k^2 < \infty$.

<caption>**Algorithm 1:** Stochastic proximal gradient with FIM preconditioning (SPG-FIM)</caption>

Require: Number of iterations $K \ge 1$; sequence of step-size $(\gamma_k)_{k \ge 1}$
1 Initialize Starting point $\theta_0 \in \mathbb{R}^d, \Delta_0$
2 for $k = 1$ to $K$ do
3 &nbsp;&nbsp;&nbsp;&nbsp;• Simulation step :
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $Z^{(k)}$ using a single step of a Hastings Metropolis procedure
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Gradient computation :
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $v_k = \frac{1}{N} \sum_{i=1}^N \nabla \log f_{\theta_k}(\mathcal{D}_i, Z_i^{(k)})$
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• FIM computation :
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Compute the stochastic approximation
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\forall i \in \{1, \dots, N\}$,
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\Delta_i^{(k)} = (1 - \gamma_k)\Delta_i^{(k-1)} + \gamma_k \nabla \log f_{\theta_k}(\mathcal{D}_i, Z_i^{(k)})$
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Compute the FIM :
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$FIM_k = \frac{1}{N} \sum_{i=1}^N \Delta_i^{(k)} (\Delta_i^{(k)})^T$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Gradient descent :
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Forward step : $\omega_{k+1} = \theta_k - \gamma_k FIM_k^{-1} v_k$
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Backward step : $\theta_{k+1} = \mathrm{Prox}_{\mathrm{pen}_{\lambda}}(\omega_{k+1})$
16 end
17 return $\hat{\theta} = \theta_K$

Note that we can differentiate between the two step-size sequences involved in either stochastic approximation or gradient descent. However, for the sake of clarity, we have noted them in the same way here.

## 3.4 The Estimation Procedure

As explained, the proximal operator (9) has a shrinking effect on the estimator after its application, meaning that the values found for $\beta$ are smaller than expected. Therefore, the estimator of $\beta$ is biased; we thus divide the inference into two steps. An exploratory one that allows us to select the support of the vector $\beta$ through a Lasso penalization estimation and a second step of inference without penalization, where we have restricted the number of covariates with respect to the selected support. We also need to select a well balance value for the regularization parameter, for that we use the BIC criterion. We thus apply these two steps for a grid of values for the regularization parameter $\lambda$ and select the model and estimation that minimizes the BIC criterion. We conduct the following inference methodology :
where $\lambda$ is a positive real called regularization parameter. Let us define for each value of $\lambda$ in a grid $\Lambda$ the penalized maximum likelihood estimator by:

$$ \hat{\theta}_{\text{pen}}^{\lambda}(\mathcal{D}) = \arg \max_{\theta \in \Theta} (\log \mathcal{L}_{\text{marg}}(\theta; \mathcal{D}) - \text{pen}_{\lambda}(\theta)), \quad (6) $$

where $\Theta$ denotes the parameter space and $\lambda$ is a positive parameter. The larger the value of $\lambda$, the more $\beta$ will be constrained to have zero components. Conversely, the smaller the value of $\lambda$, the freer the components of $\beta$ will be. It is customary to determine the value of $\lambda$ using cross-validation (Tibshirani, 1996), in particular for prediction objective. But here, without any predictive context, we choose to fix an optimal value for the regularization parameter based on a model criterion as detailed in the next section.

### 3.3 Regularization path procedure

To select only the most explanatory covariates, it's important to choose a well-balanced value for the regularization parameter $\lambda$. We choose an optimal parameter by minimizing the extended Bayesian Information Criterion (eBIC) (see Chen and Chen, 2008). Following the study of Delattre, Lavielle, and Poursat, 2014, we penalize the number of degrees of freedom by the logarithm of the total number of observations $N \times J$. We consider a grid $\Lambda$ of values for the regularization parameter $\lambda$. We conduct then the following inference methodology :

i) For all $\lambda \in \Lambda$ repeat the following steps:

* calculate $\hat{\theta}_{\text{pen}}^{\lambda}(\mathcal{D}) = \arg \max_{\theta \in \Theta} (\log \mathcal{L}_{\text{marg}}(\theta; \mathcal{D}) - \text{pen}_{\lambda}(\theta)),$

* deduce the associated support $\hat{S}_{\lambda} = \{j \in \{1, \dots, p\}, \hat{\beta}_{\lambda,j}^{\text{PEN}} \neq 0\}$

* calculate $\hat{\theta}_{\lambda}^{\text{MLE}}(\mathcal{D}) = \arg \max_{\theta \in \Theta} \{\log g_{\lambda}(\theta; \mathcal{D})\}$, where $g_{\lambda}$ is the marginal likelihood in the model reduced to the support $\hat{S}_{\lambda}$.

* calculate the eBIC criterion:

$$ \text{eBIC}(\lambda) = -2 \log g_{\lambda}(\hat{\theta}_{\lambda}^{\text{MLE}}; Y) + |\hat{S}_{\lambda}| \log(NJ) + 2 \log \left( \frac{p}{|\hat{S}_{\lambda}|} \right) \quad (7) $$

ii) Select the parameter $\lambda$ that minimizes the eBIC:

$$ \hat{\lambda} = \arg \min_{\lambda \in \Lambda} \text{eBIC}(\lambda), $$

and consider the final estimator defined by $\hat{\theta}_{\lambda}^{\text{MLE}}$.

## 4 Numerical methods for inference procedure

To implement our variable selection and re-estimation procedures in practice, we must be able to compute both the maximum likelihood estimate in a reduced model and the penalized estimate across a range of regularization parameter values, with both efficiency and accuracy. In both cases, we face the challenge of integrating over latent variables, which typically precludes closed-form expressions for the optimization criteria. As a result, we turn to numerical methods to solve these maximization problems. Specifically, we propose using an adaptive stochastic gradient algorithm for the first objective and a proximal extended variant for the second.
* Run the SPG-FIM procedure (algorithm 1) in order to compute

$$ \hat{\theta}_{\text{Lasso}}^{\lambda}(\mathbf{D}) = \arg \max_{\theta \in \Theta} \{\log \mathcal{L}_{\text{marg}}(\theta; \mathbf{D}) - \text{pen}_{\lambda}(\theta)\}, $$

for different values of $\lambda$ on a fixed grid.

* Choose the reduced support of $\beta$ according to the estimate $\hat{\theta}_{\text{Lasso}}^{\lambda}(\mathbf{D})$. We got a reduced model. Compute $\hat{\theta}_{\text{MLE}}^{\lambda}(\mathbf{D})$ the maximum likelihood estimate in this reduced model without the penalty term and therefore without bias with the SG-FIM algorithm (SPG-FIM without the *Backward* step).

* We obtain a collection of estimates for each possible value of $\lambda$. For each regularization parameter, we compute the BIC criterion (see [23]):

$$ \text{BIC}(\lambda) = -2 \log(\mathcal{L}_{\text{marg}}(\hat{\theta}_{\text{MLE}}^{\lambda}(\mathbf{D}); \mathbf{D})) + k \log(N \times (J + 1)), $$

where $k$ is the number of non-zeros components in $\beta$.

Note that the quantity $\mathcal{L}_{\text{marg}}(\theta; \mathbf{D}) = \int_{\mathbf{Z}} f_{\theta}(\mathbf{D}, \mathbf{Z})d\mathbf{Z}$ is computed by approximating the integral using a Monte Carlo procedure.

* Finally, we keep the estimate $\hat{\theta}_{\text{MLE}}^{\lambda}(\mathbf{D})$ and the regularization parameter $\hat{\lambda}$ that minimizes the BIC, as :

$$ \hat{\lambda} = \arg \min_{\lambda} \text{BIC}(\lambda). $$

So, with this method, we can finally estimate the relevant parameters from the joint modeling of survival and longitudinal data.

# 4 Simulation study

In this section, we propose to study the performance of the procedure we have just presented. First, we want to numerically show the consistency of the estimator (7) by studying four scenarios, varying the number $N$ of individuals observed. Then, we will demonstrate the method’s ability to select variables by varying the number $p$ of covariates in the initial model. This simulation study is conducted using joint modeling for longitudinal logistic data. In each scenario, we generate an independent data set, and each time, we fit the model (4) using the routine described in subsection 3.4.

To compare the results of the different scenarios, we look at two metrics: the root mean square error (RMSE) to measure the estimation quality of the method and the sensitivity, specificity, and accuracy to study the selection capacity of the method. Sensitivity is the proportion of true positives correctly identified. Specificity is the proportion of the true negatives correctly identified. Finally, accuracy is the proportion of true results, either true positive or true negative.

## 4.1 Logistic Joint-Model Specification

We consider the joint model defined in section 2 with the classical logistic function for the nonlinear mixed-effects model, defined by :

$$ m: t \mapsto \frac{\mathbf{Z}_1}{1 + \exp\left(\frac{\mathbf{Z}_2 - t}{\tau}\right)}, \qquad (10) $$

where $\mathbf{Z}_1$ represents the asymptotical maximum value of the curve, $\mathbf{Z}_2$ represents the value of the sigmoid’s midpoint, and $\tau$ represents the logistic growth rate. We model for each individual $i$ the corresponding individual parameter $\mathbf{Z}_i \in \mathbb{R}^2$ through a Gaussian random variable with expectation $\mu \in \mathbb{R}^2$ and diagonal variance-covariance matrix $\Gamma = \text{diag}(\gamma_1^2, \gamma_2^2)$. Regarding the survival model, we consider a fixed Weibull
## 4.1 Adaptive Stochastic Gradient Descent Algorithm

Adaptive algorithms such as AdaGrad (Duchi, Hazan, and Singer, 2011), RMSProp (Tieleman, 2012), and Adam (Kingma and Ba, 2017) have demonstrated their effectiveness. These methods aim to automatically rescale each direction of the gradient descent by using an adaptive learning rate, thereby adjusting the scale of individual gradient components and ensuring a more balanced evolution of the optimization process.

In our work, we choose to leverage the advantages of the AdaGrad algorithm. Its primary strength lies in the adaptivity of its learning rate, which is crucial since the convergence behavior of stochastic gradient algorithms is highly sensitive to this choice. Additionally, AdaGrad's component-wise independence in adjusting the learning rate proves particularly beneficial when handling the non-differentiable penalty term involved in the criterion discussed in the next section.

The algorithm proceeds in two main steps. First, we simulate a realization of the latent variables, either by directly sampling from their posterior distribution using the current parameter values or by performing a step of a Metropolis-Hastings sampler (Marin and Robert, 2007). In the second step, the parameter estimates are updated via adaptive preconditioned gradient descent. More precisely, given the sampled latent variables, we compute the gradient of the complete log-likelihood and, for each parameter component, cumulate the squared gradients over iterations. These cumulative values are then used to define a preconditioning diagonal matrix, with entries inversely proportional to the square root of the cumulated gradients.

## 4.2 Adaptive Weighted Proximal Stochastic Gradient Descent Algorithm

Due to the non-differentiability of the considered $l_1$-penalty term, we use a proximal algorithm as presented by Achab, 2017 and Fort, Ollier, and Samson, 2019. Thus, we add a proximal step in the adaptive stochastic gradient algorithm to calculate the penalized estimate defined in Equation (6).

The algorithm is divided into three steps; a realization of the latent variables is sampled with a first step called *Simulation*, which uses a Metropolis-Hastings sampler (see Marin and Robert, 2007). The second step is the classical gradient descent on the gradient of the logarithm of the complete likelihood, the *Forward* step. The last step, called *Backward*, deals with the penalty term. We apply the classical proximal operator (see Moreau, 1962; Rockafellar, 1976), defined below:

$$ \mathrm{Prox}_{\mathrm{pen}_{\lambda}}(\beta) = \arg\min_{\beta' \in \mathbb{R}^p} \left( \mathrm{pen}_{\lambda}(\beta') + \frac{1}{2} \|\beta - \beta'\|_2^2 \right). \quad (8) $$

Considering the $l_1$-penalty, the proximal operator has an explicit form given for all $i \in \{1, \dots, p\}$:

$$ (\mathrm{Prox}_{\lambda} \cdot |(\beta))_i = \begin{cases} 0 & \text{if } |\beta_i| < \lambda \\ \beta_i - \lambda & \text{if } \beta_i \ge \lambda \\ \beta_i + \lambda & \text{if } \beta_i \le -\lambda \end{cases} \quad (9) $$

The *Backward* step corresponds to the application of the proximal operator on the result of the *Forward* step.

# 5 Numerical experiments

In this section, we study the performance of the proposed procedure called NLJM-LASSO. We assess its performance in different scenarios. We first study the impact of the censoring rate on the variable selection and estimation of the parameters when individual parameters of the mixed model are independent. Then, we will highlight the difficulty appearing when the individual parameters of
baseline defined as $h_{a,b}(t) = ba^{-b}t^{b-1}$, where $a = 800$ and $b = 10$ are fixed (i.e. not estimated) in the simulation study. Moreover we assume that the time of interest $T^*$ can be related to the longitudinal dynamic using the following linking function: $f(\alpha, \mathcal{M}(t; \mathbf{Z}_i)) = \alpha m(t, \mathbf{Z}_i)$.

We focus in this simulation study on the variable selection procedure and on the inference of the parameters of the mixed-effects model as well as $\alpha$, the multiplicative parameter of the Cox model.

### 4.1.1 Simulation setup

We generated 100 data sets according to the joint model presented above. For each different value of $p$, we choose the vector $\beta$ such that the first four components are equal to $(-3, -2, 2, 3)$ and the rest are equal to zero. Additionally, we generate the matrix of covariates $U$ with $N$ rows and $p$ columns, following a uniform distribution $U_{i,l} \sim U([-1, 1])$; $\forall i \in 1, ..., N, l \in 1, ..., p$. Afterward, we renormalize all the matrix columns to obtain centered data. All the parameter values are detailed in the table 4.1.1 below.

<table><caption>Table 1: Parameter values used in simulation</caption><thead><tr><th>Parameter</th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>τ</th><th>γ<sub>1</sub><sup>2</sup></th><th>γ<sub>2</sub><sup>2</sup></th><th>σ<sup>2</sup></th><th>α</th></tr></thead><tbody><tr><td>True value</td><td>200</td><td>500</td><td>150</td><td>40</td><td>100</td><td>10</td><td>0.05</td></tr><tr><td>Parameter</td><td>β<sub>1</sub></td><td>β<sub>2</sub></td><td>β<sub>3</sub></td><td>β<sub>4</sub></td><td>β<sub>5</sub></td><td>...</td><td>β<sub>p</sub></td></tr><tr><td>True value</td><td>-3</td><td>-2</td><td>2</td><td>3</td><td>0</td><td>...</td><td>0</td></tr></tbody></table>

The initial parameters $\theta_0$ were randomly drawn for each independent run. To make the algorithm numerically more robust, we add a regularization to the Fisher information matrix when estimating it. The Fisher information is thus calculated as follows: $FIM_k = (1 - \lambda_k)Id + \lambda_k \frac{1}{N} \sum_{i=1}^{N} \Delta_i^{(k)} (\Delta_i^{(k)})^T$. As explained before, we can set up several step-size sequences for stochastic gradient approximation, gradient descent, and fisher regularization. We define them as follows $\lambda_{sto}, \lambda_{grad}$ and $\lambda_{fim}$. The gradient descent step size is chosen to be first increasing, then stable for a given time interval, and then decreasing. This provides a warm-up period to stabilize the Gibbs sampler. During this warm-up phase, step sizes for the stochastic approximation and the regularization of the Fisher information were also chosen to increase and then become stable until the end of the procedure. We choose the following sequences of steps according to $k_{preheating}, k_{heating}$ the start and end iterations of the preheating and heating phases respectively :

$$ \gamma_{sto,k} = \begin{cases} \gamma_0 \left(1 - \frac{k}{k_{preheating}}\right) & \text{if } k \le k_{preheating} \\ 1 & \text{else} \end{cases} $$ $$ \gamma_{grad,k} = \begin{cases} \gamma_0 \left(1 - \frac{k}{k_{preheating}}\right) & \text{if } k \le k_{preheating} \\ 1 & \text{else if } k \le k_{heating} \\ (k - k_{heating})^{-0.65} & \text{else} \end{cases} \quad (11) $$ $$ \gamma_{fim,k} = \begin{cases} \gamma_0 \left(1 - \frac{k}{k_{preheating}}\right) & \text{if } k \le k_{preheating} \\ 0.9 & \text{else} \end{cases} $$

In this simulation study, we use the following configuration: $k_{preheating} = 1000$, $k_{heating} = 1400$ and $\gamma_0 = 10^{-9}$.

Table 2 summarizes the root mean square errors (RMSE) obtained when the number of covariates $p$ is equal to 500, the number of observations $J$ is fixed at 15, and the number of individuals increases $N \in \{200, 400\}$. We aim to show the numerical consistency of the proposed estimator (7). A very small error is obtained for a large proportion of the parameters. The results improve as the number of individuals increases. The variance of random effects is difficult to estimate, particularly that associated with the position of the midpoint in the logistic function $\gamma_2^2$. It can be seen that the parameters related to the explanatory covariates (i.e. $\beta_1, \beta_2, \beta_3, \beta_4$) are more sensitive than the other parameters to the number of individuals. Excluding the variances, the algorithm allows us to obtain estimates below an error of 5% for $N = 400$.
<caption>**Algorithm 1:** Adaptive Weighted Proximal Stochastic Gradient Algorithm</caption>

Require: $\gamma_0 > 0$, $\epsilon > 0$, $K_{max} \in \mathbb{N}$
1 Initialize Initialize starting point $\theta_0 \in \mathbb{R}^d$; adaptive stepsize vector $s_0 = 0$; iteration $k = 0$.
2 while $(\theta_k) \text{ not converged or } k \le K_{max}$ do
3 $k \leftarrow k + 1$
4 • Simulation step:
5 Draw $Z^{(k)}$ from the posterior distribution $p(\cdot|\theta_k)$ or using a single step of a Metropolis Hastings algorithm having as stationary distribution the posterior.
6 • Gradient computation: Compute $v_k = \frac{1}{N} \sum_{i=1}^N \nabla \log f_{\theta_k}(\mathbf{D}_i, \mathbf{Z}_i^{(k)})$
7 • Adaptive stepsize: $s_{k,l}^2 \leftarrow s_{k-1,l}^2 + v_{k,l}^2$ for $l \in \{1, \dots, p\}$
8 • Update parameters:
9 • Forward step: $\omega_{k,l} \leftarrow \theta_{k-1,l} + \gamma_0 v_{k,l} / (s_{k,l} + \epsilon)$ for $l \in \{1, \dots, p\}$
10 • Backward step: $\theta_{k,l} \leftarrow \text{Prox}_{(P_{k,l}+\epsilon)/\gamma_0, \lambda|\cdot|}(\omega_{k,l})$ for $l \in \{1, \dots, p\}$
11 end
12 Return $\hat{\theta} = \theta_k$

the model are correlated. In each scenario, we generate $n_{\text{run}}$ independent data sets and for each run we fit the considered model using the routine described in Section 3.3 and get an estimate $\hat{\theta}_i^{\text{MLE}}$ for each $i \in \{1, \dots, n_{\text{run}}\}$. To compare the results of the different scenarios, we evaluate the variable selection capacity of the method by computing the sensitivity, specificity and accuracy. Sensitivity is the proportion of true positives (TP) correctly identified. Specificity is the proportion of the true negatives (TN) correctly identified. Finally, accuracy is the proportion of true results, either true positive or true negative. We abbreviated false negatives and false positives respectively by FN and FP.

$$ \text{Se} = \frac{TP}{TP + FN} \quad \text{Sp} = \frac{TN}{TN + FP} \quad \text{Ac} = \frac{TP + TN}{P + N} $$

We look also at the mean square error (mse) and the relative root mean square errors (RRMSE) to measure the estimation quality of the method defined as:

$$ \text{RRMSE}(\hat{\theta}_k^{\text{MLE}}) = \sqrt{\frac{1}{n_{\text{run}}} \sum_{i=1}^{n_{\text{run}}} \frac{(\hat{\theta}_{i,k}^{\text{MLE}} - \theta_k^*)^2}{\theta_k^{*2}}} \quad \text{mse}(\hat{\theta}_k^{\text{MLE}}) = \frac{1}{n_{\text{run}}} \sum_{i=1}^{n_{\text{run}}} (\hat{\theta}_{i,k}^{\text{MLE}} - \theta_k^*)^2 $$

where $\theta^*$ denotes the true value of the parameters.

## 5.1 Variable selection in high dimension and re-estimation analysis in a non linear joint model with independent random effects from censored survival data

We consider the joint causal model defined in Section 2 with the classical logistic growth function for the non linear mixed-effect model defined by:

$$ Y_{i,j} = \frac{Z_{i1}}{1 + \exp(-(t_{i,j} - Z_{i2})/\tau)} + \varepsilon_{i,j} \quad (10) $$

where the quantity $Z_{i1}$ represents the asymptotic value of the curve, $Z_{i2}$ the value of the sigmoid midpoint, and $\tau$ the logistic growth rate. The term $\varepsilon_{i,j}$ represents an additive noise term. We
<table><caption>Table 2: Root mean squared error (RMSE) and Relative-RMSE (RMSE divided by the simulated true value), for a total of 100 repetitions. The simulated values are $\mu_1 = 200$, $\mu_2 = 500$, $\tau = 150$, $\gamma_1^2 = 40$, $\gamma_2^2 = 100$, $\sigma^2 = 10$, $\alpha = 0.05$, $\beta = [-3, -2, 2, 3]$</caption><thead><tr><th rowspan="3"></th><th colspan="4">P = 500</th></tr><tr><th colspan="2">N = 200</th><th colspan="2">N = 400</th></tr><tr><th>RMSE</th><th>RRMSE(%)</th><th>RMSE</th><th>RRMSE(%)</th></tr></thead><tbody><tr><th>μ<sub>1</sub></th><td>0.467</td><td>0.233</td><td>0.303</td><td>0.151</td></tr><tr><th>μ<sub>2</sub></th><td>0.793</td><td>0.159</td><td>0.623</td><td>0.125</td></tr><tr><th>τ</th><td>0.371</td><td>0.248</td><td>0.276</td><td>0.184</td></tr><tr><th>γ<sub>1</sub><sup>2</sup></th><td>4.596</td><td>11.489</td><td>3.105</td><td>7.763</td></tr><tr><th>γ<sub>2</sub><sup>2</sup></th><td>15.501</td><td>15.501</td><td>9.542</td><td>9.542</td></tr><tr><th>σ<sup>2</sup></th><td>0.259</td><td>2.587</td><td>0.209</td><td>2.094</td></tr><tr><th>α</th><td>0.001</td><td>1.390</td><td>0.001</td><td>1.060</td></tr><tr><th>β<sub>1</sub></th><td>0.109</td><td>3.621</td><td>0.088</td><td>2.928</td></tr><tr><th>β<sub>2</sub></th><td>0.119</td><td>5.949</td><td>0.084</td><td>4.195</td></tr><tr><th>β<sub>3</sub></th><td>0.143</td><td>7.170</td><td>0.097</td><td>4.862</td></tr><tr><th>β<sub>4</sub></th><td>0.131</td><td>4.372</td><td>0.094</td><td>3.144</td></tr></tbody></table>

In addition to the estimation, we also want to demonstrate our method's ability to perform variable selection. To do this, we set the number of individual parameters to $N = 200$ and $J = 15$ and study the cases where $p \in \{200, 300, 500\}$. Table 3 shows the results for these scenarios.

<table><caption>Table 3: Average of sensitivity (Se), specificity (Sp), accuracy (Ac), and error using $L_1$-norm and $L_2$-norm over 100 repetitions, with their associated variance in brackets.</caption><thead><tr><th rowspan="2">P</th><th colspan="5">N = 200</th></tr><tr><th>Se</th><th>Sp</th><th>Ac</th><th>||·||<sub>1</sub></th><th>||·||<sub>2</sub></th></tr></thead><tbody><tr><td>200</td><td>1.0<br>(0.0)</td><td>0.999<br>(1.62e-05)</td><td>0.999<br>(1.57e-05)</td><td>0.478<br>(0.123)</td><td>0.268<br>(0.0255)</td></tr><tr><td>300</td><td>1.0<br>(0.0)</td><td>0.999<br>(8.46e-06)</td><td>0.999<br>(8.29e-06)</td><td>0.481<br>(0.119)</td><td>0.273<br>(0.0251)</td></tr><tr><td>500</td><td>1.0<br>(0.0)</td><td>0.999<br>(2.23e-06)</td><td>0.999<br>(2.2e-06)</td><td>0.548<br>(0.131)</td><td>0.307<br>(0.0305)</td></tr></tbody></table>

We can see that sensitivity, i.e., the proportion of true positives, is 100% with zero variance, i.e., over the 100 repetitions each time the four explanatory covariates were selected. The proportion of false positives, or specificity, is almost equally 100%. For the three values of $p$ tested in less than 15% of cases, one to four additional variables are selected at worst. The error selection error for the $L_1$ or $L_2$-norms is constant for the scenarios we tested.

As an illustration, Figure 1 displays the estimated parameter as a function of iterations during the execution of the SPG-FIM algorithm. In addition to the algorithm's convergence to the values used in the simulation, we can also observe the preheating and heating phases before stabilizing the algorithm.
model for each individual i the corresponding individual parameter $Z_i \in \mathbb{R}^2$ through a Gaussian random variable with expectation $\mu = (\mu_1, \mu_2) \in \mathbb{R}^2$ and diagonal covariance matrix $\Omega = \text{diag}(\omega_1^2, \omega_2^2)$ with $\omega_1^2, \omega_2^2 \in (\mathbb{R}_+^*)$. We assume that the noise term follows a centered Gaussian distribution with unknown variance $\sigma^2 \in \mathbb{R}_+^*$. Regarding the survival model, we consider a fixed unknown constant value for the baseline $h_0(t) = h_0, \forall t \ge 0$. Moreover we assume that the time of interest $\mathbf{T}^*$ can be related to the longitudinal dynamic using the following link function: $f(\alpha, \mathcal{M}(t; Z_i)) = \alpha m(t, Z_i)$ with

$$m(t, Z_i) = \frac{Z_{i1}}{1 + \exp(-(t_{i,j} - Z_{i2})/\tau)}, \quad \forall t \ge 0 \qquad (11)$$

Thus the hazard rate writes: $\forall i \in \{1, ..., N\}, \forall j \in \{1, ..., J\}$:

$$h(t|\mathcal{M}(t, Z_i), X_i) = h_0 \exp(\beta^T X_i + \alpha m(t, Z_i)), \quad \forall t \ge 0 \qquad (12)$$

The whole joint model parameter are denoted by $\theta = (h_0, \beta, \alpha, \mu_1, \mu_2, \tau, \omega_1^2, \omega_2^2, \sigma^2)$. We generated independently 100 data sets with the following parameter's values: $h_0 = 0.05, \mu_1 = 4, \mu_2 = 15, \tau = 4, \omega_1^2 = 0.1, \omega_2^2 = 1, \sigma^2 = 0.05$, and $\alpha = 1.5$. Each vector of covariates $X_i$ is generated with p components, where each components is drawn independently from a binomial distribution with a success probability of 0.2, $X_{i,k} \sim \mathcal{B}(0.2), \forall i \in \{1..., N\}, k \in \{1..., p\}$. The first four components of the vector $\beta$ are set to $(-2, -1, 2, 2)$, while the remaining components are set to zero. We simulate censoring of the event $T_i^*$ by generating a random variable $C_i$ from a uniform distribution on the interval $[0, T_i^*]$ and a censoring indicator $\delta_i$ according to an independent Bernoulli distribution with a probability $\rho$ of being observed equal to the desired observation percentage. Each individual's survival time has then a probability of being censored equal to $1 - \rho$. The event time $T_i$ is defined by $T_i = C_i \mathbb{1}_{\delta_i=0} + T_i^* \mathbb{1}_{\delta_i=1}$. The number of individuals is fixed at $N = 200$ and the number of observations per individual is fixed at $J = 15$. The number of covariates is fixed at $p = 500$. In this simulation study, we focus on the variable selection procedure as well as on the parameter estimation in the reduced model once the covariate support has been selected in a context of censored data. Thus we consider three scenarios with different percentages of censored data: without censoring 0%, with moderate censoring 40%, and heavy censoring 70%.

The results of the selection procedure are presented in Table 1. The selection score are given, as well as the mean square error (mse). We observe that the censoring has a strong impact on the selection procedure. The sensitivity, which is the proportion of relevant covariates selected, decreases strongly as the percentage of censoring increases. Figure 1 presents the proportion of correctly selected non-zero components for the first four components of $\beta$. We observe that these proportions decrease strongly as the censoring rate increases. In all scenario the procedure struggles to select non-zero components of $\beta$ when the censoring is high. The component $\beta_2$ seems to be more difficult to select than the others. This can be explained by its relative lower amplitude.

<table><caption>Table 1: Average sensitivity (Se), specificity (Sp), accuracy (Ac), and mean square error (mse) using the procedure NLJM-LASSO over 100 repetitions for the logistic joint model.</caption><thead><tr><th></th><th>0%</th><th>40%</th><th>70%</th></tr></thead><tbody><tr><td>Ac</td><td>0.999</td><td>0.997</td><td>0.994</td></tr><tr><td>Se</td><td>0.940</td><td>0.782</td><td>0.230</td></tr><tr><td>Sp</td><td>1.000</td><td>0.999</td><td>1.000</td></tr><tr><td>mse</td><td>0.001</td><td>0.005</td><td>0.019</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Parameter estimates across the SPG-FIM iterations. Dotted line: simulated value, red vertical line: end of the pre-heating, green vertical line: end of the heating</figcaption></figure>

# 5 Conclusion and perspectives

In this work, we jointly addressed variable selection and parameter estimation in a joint model for survival data and longitudinal outcomes. We connected a survival model with a nonlinear mixed-effects model through a linking function and covariates eventually in high dimension. To estimate the unknown parameters of this joint model, we conduct a preconditioned proximal stochastic gradient to deal with the latent variables and the Lasso penalty. Our methodology has been thoroughly evaluated on simulated data to demonstrate its performance.

One interesting perspective of this work consists of addressing the prediction task carefully. It would be interesting to set up a method for predicting survival time from some start of longitudinal data observation. Previous studies, such as [8] and [13]), have explored this task within the context of joint models.

Furthermore, considering the partial likelihood for the survival part rather than the complete likelihood would be an intriguing avenue to explore. This approach would eliminate the need to make assumptions about the form of the baseline risk, which is especially relevant when dealing with the Cox model. Indeed, when considering the Cox model, the partial likelihood allows for estimating the regression parameter without knowing the baseline function. It would be interesting to see how our algorithm performs with this partial likelihood instead of the complete one.

# Funding and Acknowledgements

This work was funded by the https://stat4plant.mathnum.inrae.fr/(Stat4Plant) project ANR-20-CE45-0012.

# References

[1] Massil Achab. *Learning from Sequences with Point Processes*. PhD thesis, Université Paris-Saclay, 2017. Issue: 2017SACLX068.

[2] Baey, Delattre, Kuhn, Leger, and Lemler. Efficient preconditioned stochastic gradient descent for estimation in latent variables models. *ICML*, 2023.
<figure><img src="image_1.png" /><figcaption>Figure 1: Proportion of correctly selected non-zero components for the logistic joint model using the NLJM-LASSO procedure over 100 repetitions with three different censoring rates.</figcaption></figure>

After the selection procedure, we re-estimate the parameters in the reduced model. The results are summarized in Table 2. The censoring rate seems to have globally a small impact on the estimation of the parameter of the mixed-effect models in this configuration. The variances of the random effects are poorly estimated whatever the censoring rate. The variance $\omega_1^2$ seems to be more affected than the other parameters. The link parameter $\alpha$ is more difficult to estimate when the censoring increases (from 10% to 16% of relative error) as well as the baseline constant $h_0$. The four non-zero components are well estimated without censoring. However these estimations are deeply deteriorated under moderate and heavy censoring, in particular for the component $\beta_2$. This phenomenon is coherent with results presented in Figure 1 and can be explained by the lower proportion of correctly selected variables for this specific component.

As an illustration of the convergence of the adaptive stochastic gradient algorithm, Figure 2 displays five trajectories of the re-estimation of survival parameters $(h_0, \alpha, \beta)$ as a function of iterations during the execution with different initializations. We observe that the algorithm converges to a limit close to the true value of the parameters. We emphasize here that the value targeted by the adaptive stochastic gradient algorithm is a critical point of the marginal likelihood, which turns out to be a maximum of the marginal likelihood. We recall that the asymptotic behavior expected under regularity assumptions as the sample size increases is indeed a convergence of the maximum likelihood estimate towards the true value of the parameter. This explains the fact that the limits are not exactly equal to the true values of the parameters, but only close to them.
[3] D. R. Cox. Regression models and life-tables. *Journal of the Royal Statistical Society. Series B (Methodological)*, 34(2):187-220, 1972. Publisher: [Royal Statistical Society, Wiley].

[4] M. Davidian and D.M. Giltinan. *Nonlinear Models for Repeated Measurement Data*. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis, 1995.

[5] Vianney Debavelaere and Stéphanie Allassonnière. On the curved exponential family in the stochastic approximation expectation maximization algorithm. *ESAIM: Probability & Statistics*, 25, 2021.

[6] Maud Delattre and Estelle Kuhn. Computing an empirical Fisher information matrix estimate in latent variable models through stochastic approximation. *Computo*, 2023.

[7] Solène Desmée, France Mentré, Christine Veyrat-Follet, Bernard Sébastien, and Jeremie Guedj. Using the saem algorithm for mechanistic joint models characterizing the relationship between nonlinear psa kinetics and survival in prostate cancer patients. *Biometrics*, 73(1):305-312, 2017.

[8] Solène Desmée, France Mentré, Christine Veyrat-Follet, Bernard Sébastien, and Jérémie Guedj. Non-linear joint models for individual dynamic prediction of risk of death using hamiltonian monte carlo: application to metastatic prostate cancer. *BMC Medical Research Methodology*, 17(1):105, 2017.

[9] Gersende Fort, Edouard Ollier, and Adeline Samson. Stochastic proximal gradient algorithms for penalized mixed models, 2017.

[10] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, PAMI-6(6):721-741, 1984.

[11] Jeremie Guedj, Rodolphe Thiébaut, and Daniel Commenges. Joint modeling of the clinical progression and of the biomarkers' dynamics using a mechanistic model. *Biometrics*, 67(1):59-66, 2011.

[12] Zangdong He, Wanzhu Tu, Sijian Wang, Haoda Fu, and Zhangsheng Yu. Simultaneous variable selection for joint models of longitudinal and survival outcomes: Variable selection in joint models. *Biometrics*, 71(1):178-187, 2015.

[13] Marion Kerioui, Julie Bertrand, René Bruno, François Mercier, Jérémie Guedj, and Solène Desmée. Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models. *British Journal of Clinical Pharmacology*, 88(4):1452-1463, 2022.

[14] Agnieszka Król, Loïc Ferrer, Jean-Pierre Pignon, Cécile Proust-Lima, Michel Ducreux, Olivier Bouché, Stefan Michiels, and Virginie Rondeau. Joint model for left-censored longitudinal data, recurrent events and terminal event: Predictive abilities of tumor burden for cancer evolution with application to the ffcd 2000-05 trial. *Biometrics*, 72(3):907-916, 2016.

[15] Jean Jacques Moreau. Fonctions convexes duales et points proximaux dans un espace hilbertien. *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 255:2897-2899, 1962.

[16] Shu Kay Ng, Thriyambakam Krishnan, and Geoffrey J McLachlan. The em algorithm. *Handbook of computational statistics: concepts and methods*, pages 139-172, 2012.

[17] José C. Pinheiro and Douglas M. Bates. *Mixed-Effects Models in S and S-PLUS*. Statistics and Computing. Springer-Verlag, 2000.

[18] Elizabeth R. Brown and Joseph G. Ibrahim. A bayesian semiparametric joint hierarchical model for longitudinal and survival data. *Biometrics*, 59(2):221-228, 2003.

[19] Dimitris Rizopoulos. Jm: An r package for the joint modelling of longitudinal and time-to-event data. *Journal of statistical software*, 35:1-33, 2010.

[20] Dimitris Rizopoulos. *Joint models for longitudinal and time-to-event data: With applications in R*. CRC press, 2012.
<table><caption>Table 2: Mean parameter estimates and average relative root mean square errors (RRMSE) over 100 repetitions for the logistic joint model where $\hat{\theta}$ is the estimates in the reduced model after variable selection using the adaptive stochastic gradient algorithm. The value $\theta^*$ denotes the true value of the parameter.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">θ*</th><th colspan="2">C = 0%</th><th colspan="2">C = 40%</th><th colspan="2">C = 70%</th></tr><tr><th>θ̂</th><th>RRMSE(%)</th><th>θ̂</th><th>RRMSE(%)</th><th>θ̂</th><th>RRMSE(%)</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>4.0</td><td>4.15</td><td>8.9</td><td>4.14</td><td>10.2</td><td>4.14</td><td>15.0</td></tr><tr><td>μ<sub>2</sub></td><td>15.0</td><td>15.25</td><td>4.5</td><td>15.23</td><td>5.3</td><td>15.21</td><td>7.2</td></tr><tr><td>ω<sup>2</sup><sub>1</sub></td><td>0.1</td><td>0.226</td><td>152.2</td><td>0.255</td><td>177.9</td><td>0.282</td><td>204.8</td></tr><tr><td>ω<sup>2</sup><sub>2</sub></td><td>1.0</td><td>0.69</td><td>39.6</td><td>0.67</td><td>44.0</td><td>0.67</td><td>44.5</td></tr><tr><td>τ</td><td>4.0</td><td>4.06</td><td>4.3</td><td>4.05</td><td>4.7</td><td>4.04</td><td>5.5</td></tr><tr><td>σ<sup>2</sup></td><td>0.0</td><td>0.050</td><td>6.8</td><td>0.050</td><td>7.3</td><td>0.049</td><td>8.1</td></tr><tr><td>α</td><td>1.5</td><td>1.49</td><td>14.3</td><td>1.59</td><td>20.0</td><td>1.29</td><td>26.9</td></tr><tr><td>h<sub>0</sub></td><td>0.05</td><td>0.051</td><td>17.5</td><td>0.036</td><td>32.0</td><td>0.024</td><td>53.3</td></tr><tr><td>β<sub>1</sub></td><td>-2.0</td><td>-2.01</td><td>15.5</td><td>-1.98</td><td>25.2</td><td>-1.19</td><td>64.2</td></tr><tr><td>β<sub>2</sub></td><td>-1.0</td><td>-0.87</td><td>48.2</td><td>-0.63</td><td>72.9</td><td>-0.04</td><td>99.0</td></tr><tr><td>β<sub>3</sub></td><td>2.0</td><td>2.00</td><td>17.3</td><td>1.64</td><td>46.7</td><td>0.33</td><td>93.1</td></tr><tr><td>β<sub>4</sub></td><td>2.0</td><td>1.95</td><td>19.7</td><td>1.69</td><td>43.3</td><td>0.31</td><td>93.7</td></tr></tbody></table>

## 5.2 Variable selection in high dimension and re-estimation analysis in a non linear joint model with censored survival data in the presence of correlated individuals parameter

In this section, we investigate the effect of the presence of a correlation between the individual parameters on the variable selection and re-estimation for the joint model defined in the previous section. For this purpose, we consider that the bivariate individual parameters $Z_i$ are drawn from a correlated Gaussian random variables with a covariance matrix $\Omega = \begin{pmatrix} \omega_1^2 & \omega_{12} \\ \omega_{12} & \omega_2^2 \end{pmatrix} \in \mathbb{S}_2^{++}$. We set $\omega_1^2 = 0.1$, $\omega_2^2 = 1$, and $\omega_{12} = 0.1$, which corresponds to a correlation coefficient of 0.31 between the two individual parameters. The other parameters are set as in the previous section. We simulate independently 100 data sets with the same censoring procedure as previously and with the same censoring rates 0%, 40% and 70%. The results of the variable selection are given in Table 3. We observe the same impact of the censoring percentage on the quality of the selection. However Figure 3 shows that the different non zero components of $\beta$ are selected less accurately, since we observe in mean a loss of 10% in the proportion of correctly selected variables compared to the previous model without correlation.

<table><caption>Table 3: Average sensitivity (Se), specificity (Sp), accuracy (Ac), and mean square error (mse) using the procedure NLJM-LASSO over 100 repetitions for the logistic joint model with correlated individuals' parameters.</caption><thead><tr><th></th><th>0%</th><th>40%</th><th>70%</th></tr></thead><tbody><tr><td>Ac</td><td>0.999</td><td>0.997</td><td>0.994</td></tr><tr><td>Se</td><td>0.944</td><td>0.741</td><td>0.247</td></tr><tr><td>Sp</td><td>1.000</td><td>0.999</td><td>1.000</td></tr><tr><td>mse</td><td>0.001</td><td>0.006</td><td>0.019</td></tr></tbody></table>
[21] Dimitris Rizopoulos. 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, 2016.

[22] R. Tyrrell Rockafellar. Monotone operators and the proximal point algorithm. *SIAM Journal on Control and Optimization*, 14(5):877–898, 1976.

[23] Gideon Schwarz. Estimating the dimension of a model. *The Annals of Statistics*, 6(2):461–464, 1978. Publisher: Institute of Mathematical Statistics.

[24] Robert Tibshirani. Regression shrinkage and selection via the lasso. *Journal of the Royal Statistical Society: Series B (Methodological)*, 58(1):267–288, 1996.

[25] Robert Tibshirani. The lasso method for variable selection in the cox model. *Statistics in Medicine*, 16(4):385–395, 1997.

[26] Michael S. Wulfsohn and Anastasios A. Tsiatis. A joint model for survival and longitudinal data measured with error. *Biometrics*, 53(1):330–339, 1997.

[27] Fengting Yi, Niansheng Tang, and Jianguo Sun. Simultaneous variable selection and estimation for joint models of longitudinal and failure time data with interval censoring. *Biometrics*, 78(1):151–164, 2022.
<figure><img src="image_2.png" /><figcaption>Figure 2: Five trajectories of the adaptive stochastic gradient algorithm with different initializations for the re-estimated survival parameters in the logistic joint model without censoring. The true value of the parameters is represented by a dashed line.</figcaption></figure>

The results of the re-estimation of the parameters after the selection procedure are presented in Table 4. The presence of correlation between the individual parameters seems to have a strong impact on the estimation of the covariance matrix $\Omega$. Estimations of other parameters seem to be less affected by the presence of a correlation and have similar estimation results as in the previous section.

# 6 Conclusion and perspectives

In this work, we jointly addressed variable selection and parameter estimation in a non linear joint model for survival data and longitudinal outcomes. We connected a survival model with a non linear mixed-effects model through a link function and covariates in high dimension. To perform the variable selection procedure in the joint model, we consider as criterion the marginal log-likelihood
<figure><img src="image_3.png" /><figcaption>Figure 3: Proportion of correctly selected non-zero components for the logistic joint model with correlated individuals’ parameters using the NLJM-LASSO procedure over 100 repetitions for three different censoring rates.</figcaption></figure>

<table><caption>Table 4: Mean parameter estimates and average relative root mean square errors (RRMSE) over 100 repetitions for the logistic joint model with correlated individual’s parameters, where $\hat{\theta}$ are the estimates in the reduced model after variable selection using the adaptive stochastic gradient algorithm. The value $\theta^*$ is the true value of the parameter.</caption><thead><tr><th rowspan="2"></th><th colspan="3">C = 0%</th><th colspan="3">C = 40%</th><th colspan="3">C = 70%</th></tr><tr><th>θ*</th><th>θ̂</th><th>RRMSE(%)</th><th>θ̂</th><th>RRMSE(%)</th><th>θ̂</th><th>RRMSE(%)</th><th>θ̂</th><th>RRMSE(%)</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>4.0</td><td>4.19</td><td>10.6</td><td>4.26</td><td>14.1</td><td>4.21</td><td>16.5</td><td></td><td></td></tr><tr><td>μ<sub>2</sub></td><td>15.0</td><td>15.261</td><td>4.9</td><td>15.380</td><td>6.5</td><td>15.26</td><td>8.1</td><td></td><td></td></tr><tr><td>ω<sub>1</sub><sup>2</sup></td><td>0.1</td><td>0.36</td><td>342.5</td><td>0.42</td><td>453.1</td><td>0.55</td><td>555.1</td><td></td><td></td></tr><tr><td>ω<sub>12</sub></td><td>0.1</td><td>0.372</td><td>388.4</td><td>0.39</td><td>459.8</td><td>0.56</td><td>599.5</td><td></td><td></td></tr><tr><td>ω<sub>2</sub><sup>2</sup></td><td>1.0</td><td>1.27</td><td>67.0</td><td>1.34</td><td>90.5</td><td>1.67</td><td>120.7</td><td></td><td></td></tr><tr><td>τ</td><td>4.0</td><td>4.05</td><td>4.0</td><td>4.07</td><td>5.1</td><td>4.04</td><td>6.4</td><td></td><td></td></tr><tr><td>σ<sup>2</sup></td><td>0.0</td><td>0.05</td><td>6.8</td><td>0.05</td><td>7.6</td><td>0.049</td><td>7.7</td><td></td><td></td></tr><tr><td>α</td><td>1.5</td><td>1.50</td><td>14.8</td><td>1.53</td><td>20.0</td><td>1.31</td><td>25.8</td><td></td><td></td></tr><tr><td>h<sub>0</sub></td><td>0.05</td><td>0.05</td><td>17.7</td><td>0.037</td><td>30.2</td><td>0.02</td><td>53.8</td><td></td><td></td></tr><tr><td>β<sub>1</sub></td><td>-2.0</td><td>-2.01</td><td>16.4</td><td>-1.90</td><td>25.1</td><td>-1.235</td><td>60.9</td><td></td><td></td></tr><tr><td>β<sub>2</sub></td><td>-1.0</td><td>-0.895</td><td>45.5</td><td>-0.63</td><td>72.7</td><td>-0.09</td><td>98.4</td><td></td><td></td></tr><tr><td>β<sub>3</sub></td><td>2.0</td><td>2.00</td><td>17.1</td><td>1.478</td><td>55.1</td><td>0.346</td><td>92.7</td><td></td><td></td></tr><tr><td>β<sub>4</sub></td><td>2.0</td><td>1.94</td><td>20.9</td><td>1.54</td><td>51.1</td><td>0.29</td><td>94.5</td><td></td><td></td></tr></tbody></table>

reguralized with a $\ell_1$-penalty term and run an adaptive weighted proximal stochastic gradient algorithm to deal simultaneously with the latent variables and the non-differentiability of the penalty. Then we rely on an eBIC criterion to choose an optimal value for the regularization parameter. We finally re-estimate the parameters in the reduced model with an adaptive stochastic gradient descent
algorithm. Our methodology has been thoroughly evaluated on simulated data to demonstrate
its performance. We evaluate the impact of censoring and of the presence of correlation between
individual parameters.

One interesting perspective of this work consists of addressing the prediction task carefully. It
would be interesting to set up a method for predicting survival time from some start of longitudinal
data observation. Previous studies, such as Desmée et al., 2017 and Kerioui et al., 2022, have explored
this task within the context of joint models without high dimensional covariates. Furthermore,
considering the partial likelihood for the survival part rather than the parametric complete likelihood
would be an interesting avenue to explore. This approach would eliminate the need to make
assumptions about the form of the baseline risk, which is especially relevant when dealing with the
Cox model. Indeed, when considering the Cox model, the partial likelihood allows for estimating
consistently the regression parameter without knowing the baseline function. It would be interesting
to adapt the proposed procedure to the partial likelihood instead of the parametric complete one.

Funding and Acknowledgements

This work was funded by the https://stat4plant.mathnum.inrae.fr/(Stat4Plant) project ANR-20-CE45-0012.

References

Achab, Massil (2017). "Learning from Sequences with Point Processes". Issue: 2017SACLX068. PhD thesis. Université Paris-Saclay.

Baey, Charlotte et al. (2023). "Efficient preconditioned stochastic gradient descent for estimation in latent variable models". In: *International Conference on Machine Learning*. PMLR, pp. 1430–1453.

Cappé, Olivier, Eric Moulines, and Tobias Rydén (2005). *Springer series in statistics*.

Chen, Jiahua and Zehua Chen (2008). "Extended Bayesian Information Criteria for Model Selection with Large Model Spaces". In: *Biometrika* 95.3, pp. 759–771. ISSN: 00063444, 14643510.

Cox, D. R. (1972). "Regression Models and Life-Tables". In: *Journal of the Royal Statistical Society. Series B (Methodological)* 34.2. Publisher: [Royal Statistical Society, Wiley], pp. 187–220. ISSN: 00359246.

Davidian, M. and D.M. Giltinan (1995). *Nonlinear Models for Repeated Measurement Data*. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis. ISBN: 978-0-412-98341-2.

Debavelaere, Vianney and Stéphanie Allassonnière (2021). “On the curved exponential family in the Stochastic Approximation Expectation Maximization Algorithm.” In: *ESAIM: Probability & Statistics* 25.

Delattre, Maud, Marc Lavielle, and Marie-Anne Poursat (2014). “A note on BIC in mixed-effects models”. In: *Electronic Journal of Statistics* 8.1, pp. 456–475. DOI: 10.1214/14-EJS890.

Desmée, Solène et al. (2017). “Using the SAEM algorithm for mechanistic joint models characterizing the relationship between nonlinear PSA kinetics and survival in prostate cancer patients”. In: *Biometrics* 73.1, pp. 305–312.

Desmée, Solène et al. (2017). “Nonlinear joint models for individual dynamic prediction of risk of death using Hamiltonian Monte Carlo: application to metastatic prostate cancer”. In: *BMC Medical Research Methodology* 17.1, p. 105. ISSN: 1471-2288. DOI: 10.1186/s12874-017-0382-9.
Duchi, John, Elad Hazan, and Yoram Singer (2011). "Adaptive Subgradient Methods for Online Learning and Stochastic Optimization". In: *Journal of Machine Learning Research* 12.61, pp. 2121–2159.

Fort, Gersende, Edouard Ollier, and Adeline Samson (2019). "Stochastic proximal-gradient algorithms for penalized mixed models". In: *Statistics and Computing* 29, pp. 231–253.

Guedj, Jeremie, Rodolphe Thiébaut, and Daniel Commenges (2011). "Joint modeling of the clinical progression and of the biomarkers' dynamics using a mechanistic model". In: *Biometrics* 67.1, pp. 59–66.

He, Zangdong et al. (2015). "Simultaneous variable selection for joint models of longitudinal and survival outcomes: Variable Selection in Joint Models". In: *Biometrics* 71.1, pp. 178–187. ISSN: 0006341X. DOI: 10.1111/biom.12221.

Kerioui, Marion et al. (2022). "Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models". In: *British Journal of Clinical Pharmacology* 88.4, pp. 1452–1463. ISSN: 0306-5251, 1365-2125. DOI: 10.1111/bcp.15200.

Kingma, Diederik P. and Jimmy Ba (2017). *Adam: A Method for Stochastic Optimization*. arXiv: 1412.6980 [cs.LG].

Król, Agnieszka et al. (2016). "Joint model for left-censored longitudinal data, recurrent events and terminal event: Predictive abilities of tumor burden for cancer evolution with application to the FFCD 2000-05 trial". In: *Biometrics* 72.3, pp. 907–916.

Marin, Jean-Michel and Christian P Robert (2007). *Bayesian core: a practical approach to computational Bayesian statistics*. Springer.

Moreau, Jean Jacques (1962). "Fonctions convexes duales et points proximaux dans un espace hilbertien". In: *Comptes rendus hebdomadaires des séances de l'Académie des sciences* 255, pp. 2897–2899.

Ng, Shu Kay, Thriyambakam Krishnan, and Geoffrey J McLachlan (2012). "The EM algorithm". In: *Handbook of computational statistics: concepts and methods*, pp. 139–172.

Pinheiro, José C. and Douglas M. Bates (2000). *Mixed-Effects Models in S and S-PLUS*. Statistics and Computing. New York: Springer-Verlag. ISBN: 978-0-387-98957-0. DOI: 10.1007/b98882.

R. Brown, Elizabeth and Joseph G. Ibrahim (2003). "A Bayesian semiparametric joint hierarchical model for longitudinal and survival data". In: *Biometrics* 59.2, pp. 221–228.

Rizopoulos, Dimitris (2010). "JM: An R package for the joint modelling of longitudinal and time-to-event data". In: *Journal of statistical software* 35, pp. 1–33.

— (2012). *Joint models for longitudinal and time-to-event data: With applications in R*. CRC press.

— (2016). "The R Package JMBayes for Fitting Joint Models for Longitudinal and Time-to-Event Data Using MCMC". In: *Journal of Statistical Software* 72.7, 1–46. DOI: 10.18637/jss.v072.i07.

Rockafellar, R. Tyrrell (1976). "Monotone Operators and the Proximal Point Algorithm". In: *SIAM Journal on Control and Optimization* 14.5, pp. 877–898. ISSN: 0363-0129, 1095-7138. DOI: 10.1137/0314056.

Tibshirani, Robert (1996). "Regression Shrinkage and Selection Via the Lasso". In: *Journal of the Royal Statistical Society: Series B (Methodological)* 58.1, pp. 267–288. ISSN: 00359246. DOI: 10.1111/j.2517-6161.1996.tb02080.x.

— (1997). "The lasso method for variable selection in the cox model". In: *Statistics in Medicine* 16.4, pp. 385–395. ISSN: 0277-6715, 1097-0258. DOI: 10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3.

Tieleman, T. (2012). *Lecture 6.5-rmsprop: Divide the Gradient by a Running Average of Its Recent Magnitude*.

Wulfsohn, Michael S. and Anastasios A. Tsiatis (1997). "A Joint Model for Survival and Longitudinal Data Measured with Error". In: *Biometrics* 53.1, pp. 330–339. ISSN: 0006341X, 15410420.
Yi, Fengting, Niansheng Tang, and Jianguo Sun (2022). “Simultaneous variable selection and estimation for joint models of longitudinal and failure time data with interval censoring”. In: Biometrics 78.1, pp. 151–164.

# Appendix

## Algorithm 2: Simulation Algorithm

Require: number of individuals N, number of longitudinal observation J, parameter values $\sigma^2, \beta$

1. **Generate :**

• Set $t_j = t_{max} + j \times \frac{t_{max}-t_{min}}{J} \quad \forall j \in \{1, ..., J\}$

• Generate longitudinal observation :

– **draw** $\epsilon \sim N(0, \sigma^2)$,

– **draw** $Z_i \sim N(\mu, \Gamma)$,

– **compute** $\mathbf{Y}_{i,j} = m(t_j; Z_i) + \varepsilon_{i,j}$.

• Generate covariates :

– **draw** $U_i^* \sim U([-1, 1])$,

– **define** $U_i = U_i^* - \bar{U}_i^*$, where $\bar{U}_i^*$ is the mean vector along the columns.

• Generate Survival times : $\forall i \in \{1, ..., N\}$

– **generate covariates** $U_i^* \sim U([-1, 1])$,

– **draw** $u_i \sim U([0, -1])$,

– **solve** : $T_i^* = F_i^{-1}(u_i)$ where $F_i(t) = 1 - S_i(t)$ where $S_i$ is the survival function
$S_i(t) = \exp(-\int_0^t h_i(t|M(t, Z_i))dt)$

• Censoring data

– set $C_i$ such that $\frac{\#\{T_i|T_i \le C_i\}}{N}$ is the right percentage of censoring

– **censor by setting** $T_i = \min(T_i^*, C_i)$

– **compute** $\delta_i = \mathbb{1}_{T_i^* \le C_i}$

– **remove all longitudinal observations after the time** $C_i$